Обратная билинейная интерполяция?

Некоторые ответы слегка неверно истолковали вашу вопрос. то есть. они предполагают, что вам дано значение неизвестной интерполированной функции, а не интерполированное положение p (x, y) внутри квадра, в котором вы хотите найти координаты (s, t). Это более простая проблема, и гарантированно найдется решение, которое будет пересекать две прямые линии через квад.

Одна из линий будет прорезать отрезки p0p1 и p2p3, а другая - через p0p2 и p1p3, похож на выровненный по оси корпус. Эти линии однозначно определяются положением p (x, y) и, очевидно, будут пересекаться в этой точке.

Рассматривая только линию, проходящую через p0p1 и p2p3, мы можем представить семейство таких линий для каждой отдельной линии. Мы выбираем s-значение, каждое из которых разрезает квад по разной ширине. Если мы зафиксируем s-значение, мы можем найти две конечные точки, установив t = 0 и t = 1.

Итак, сначала предположим, что строка имеет вид: y = a0 * x + b0

Тогда мы знаем две конечные точки этой линии, если мы фиксируем данное значение s. Это:

(1-s) p0 + (s) p1

(1-s) p2 + (s) p3

Учитывая эти две конечные точки, мы можем определить семейство линий, подключив их в уравнение для прямой и решения для a0 и b0 как функции s . Установка s-значения дает формулу для конкретной строки. Все, что нам сейчас нужно, это выяснить, какая линия в этом семействе попадает в нашу точку p (x, y). Просто вставив координаты p (x, y) в нашу линейную формулу, мы можем затем найти целевое значение s.

Соответствующий подход можно сделать и для нахождения t.

Since you're working in 2D, your bilerp function is really 2 equations, 1 for x and 1 for y. They can be rewritten in the form:

x = t * s * A.x + t * B.x + s * C.x + D.x
y = t * s * A.y + t * B.y + s * C.y + D.y

Where:

A = p3 - p2 - p1 + p0
B = p2 - p0
C = p1 - p0
D = p0

Rewrite the first equation to get t in terms of s, substitute into the second, and solve for s.

Если все, что у вас есть, это одно значение p, такое, что p находится между минимальным и максимальным значениями в четырех углах квадрата, то нет, в общем случае найти невозможно ЕДИНСТВЕННОЕ решение (s, t), такое, что билинейный интерполант даст вам это значение.

В общем, внутри квадрата будет бесконечное число решений (s, t). Они будут лежать по изогнутой (гиперболической) траектории через квадрат.

Если ваша функция является векторной, значит, у вас есть два известных значения в некоторой неизвестной точке в квадрате? Учитывая известные значения двух параметров в каждом углу квадрата, решение МОЖЕТ существовать, но в этом нет уверенности. Помните, что мы можем рассматривать это как две отдельные, независимые проблемы. Решение каждого из них будет лежать вдоль гиперболической линии контура через квадрат. Если пара контуров пересекается внутри квадрата, то решение существует. Если они не пересекаются, значит, решения не существует.

Вы также спрашиваете, существует ли простая формула для решения проблемы. Извините, но не совсем то, что я вижу. Как я уже сказал, кривые являются гиперболическими.

Одним из решений является переключение на другой метод интерполяции. Поэтому вместо билинейного делите квадрат на пару треугольников. Внутри каждого треугольника теперь мы можем использовать действительно линейную интерполяцию. Итак, теперь мы можем решить линейную систему из 2 уравнений с 2 ​​неизвестными в каждом треугольнике. В каждом треугольнике может быть одно решение, за исключением редкого вырожденного случая, когда соответствующие кусочно-линейные контурные линии оказываются со-падающими.

Вы также спрашиваете, существует ли простая формула для решения проблемы. Извините, но не совсем то, что я вижу. Как я уже сказал, кривые являются гиперболическими.

Одним из решений является переключение на другой метод интерполяции. Поэтому вместо билинейного делите квадрат на пару треугольников. Внутри каждого треугольника теперь мы можем использовать действительно линейную интерполяцию. Итак, теперь мы можем решить линейную систему из 2 уравнений с 2 ​​неизвестными в каждом треугольнике. В каждом треугольнике может быть одно решение, за исключением редкого вырожденного случая, когда соответствующие кусочно-линейные контурные линии оказываются со-падающими.

Вы также спрашиваете, существует ли простая формула для решения проблемы. Извините, но не совсем то, что я вижу. Как я уже сказал, кривые являются гиперболическими.

Одним из решений является переключение на другой метод интерполяции. Поэтому вместо билинейного делите квадрат на пару треугольников. Внутри каждого треугольника теперь мы можем использовать действительно линейную интерполяцию. Итак, теперь мы можем решить линейную систему из 2 уравнений с 2 ​​неизвестными в каждом треугольнике. В каждом треугольнике может быть одно решение, за исключением редкого вырожденного случая, когда соответствующие кусочно-линейные контурные линии оказываются со-падающими.

Теперь мы можем использовать действительно линейную интерполяцию. Итак, теперь мы можем решить линейную систему из 2 уравнений с 2 ​​неизвестными в каждом треугольнике. В каждом треугольнике может быть одно решение, за исключением редкого вырожденного случая, когда соответствующие кусочно-линейные контурные линии оказываются со-падающими.

Теперь мы можем использовать действительно линейную интерполяцию. Итак, теперь мы можем решить линейную систему из 2 уравнений с 2 ​​неизвестными в каждом треугольнике. В каждом треугольнике может быть одно решение, за исключением редкого вырожденного случая, когда соответствующие кусочно-линейные контурные линии оказываются со-падающими.

Вот моя реализация решения Naaff, как вики сообщества. Еще раз спасибо.

Это реализация C, но она должна работать на c ++. Он включает в себя функцию нечеткого тестирования.

---

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

int equals( double a, double b, double tolerance )
{
    return ( a == b ) ||
      ( ( a <= ( b + tolerance ) ) &&
        ( a >= ( b - tolerance ) ) );
}

double cross2( double x0, double y0, double x1, double y1 )
{
    return x0*y1 - y0*x1;
}

int in_range( double val, double range_min, double range_max, double tol )
{
    return ((val+tol) >= range_min) && ((val-tol) <= range_max);
}

/* Returns number of solutions found.  If there is one valid solution, it will be put in s and t */
int inverseBilerp( double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3, double x, double y, double* sout, double* tout, double* s2out, double* t2out )
{
    int t_valid, t2_valid;

    double a  = cross2( x0-x, y0-y, x0-x2, y0-y2 );
    double b1 = cross2( x0-x, y0-y, x1-x3, y1-y3 );
    double b2 = cross2( x1-x, y1-y, x0-x2, y0-y2 );
    double c  = cross2( x1-x, y1-y, x1-x3, y1-y3 );
    double b  = 0.5 * (b1 + b2);

    double s, s2, t, t2;

    double am2bpc = a-2*b+c;
    /* this is how many valid s values we have */
    int num_valid_s = 0;

    if ( equals( am2bpc, 0, 1e-10 ) )
    {
        if ( equals( a-c, 0, 1e-10 ) )
        {
            /* Looks like the input is a line */
            /* You could set s=0.5 and solve for t if you wanted to */
            return 0;
        }
        s = a / (a-c);
        if ( in_range( s, 0, 1, 1e-10 ) )
            num_valid_s = 1;
    }
    else
    {
        double sqrtbsqmac = sqrt( b*b - a*c );
        s  = ((a-b) - sqrtbsqmac) / am2bpc;
        s2 = ((a-b) + sqrtbsqmac) / am2bpc;
        num_valid_s = 0;
        if ( in_range( s, 0, 1, 1e-10 ) )
        {
            num_valid_s++;
            if ( in_range( s2, 0, 1, 1e-10 ) )
                num_valid_s++;
        }
        else
        {
            if ( in_range( s2, 0, 1, 1e-10 ) )
            {
                num_valid_s++;
                s = s2;
            }
        }
    }

    if ( num_valid_s == 0 )
        return 0;

    t_valid = 0;
    if ( num_valid_s >= 1 )
    {
        double tdenom_x = (1-s)*(x0-x2) + s*(x1-x3);
        double tdenom_y = (1-s)*(y0-y2) + s*(y1-y3);
        t_valid = 1;
        if ( equals( tdenom_x, 0, 1e-10 ) && equals( tdenom_y, 0, 1e-10 ) )
        {
            t_valid = 0;
        }
        else
        {
            /* Choose the more robust denominator */
            if ( fabs( tdenom_x ) > fabs( tdenom_y ) )
            {
                t = ( (1-s)*(x0-x) + s*(x1-x) ) / ( tdenom_x );
            }
            else
            {
                t = ( (1-s)*(y0-y) + s*(y1-y) ) / ( tdenom_y );
            }
            if ( !in_range( t, 0, 1, 1e-10 ) )
                t_valid = 0;
        }
    }

    /* Same thing for s2 and t2 */
    t2_valid = 0;
    if ( num_valid_s == 2 )
    {
        double tdenom_x = (1-s2)*(x0-x2) + s2*(x1-x3);
        double tdenom_y = (1-s2)*(y0-y2) + s2*(y1-y3);
        t2_valid = 1;
        if ( equals( tdenom_x, 0, 1e-10 ) && equals( tdenom_y, 0, 1e-10 ) )
        {
            t2_valid = 0;
        }
        else
        {
            /* Choose the more robust denominator */
            if ( fabs( tdenom_x ) > fabs( tdenom_y ) )
            {
                t2 = ( (1-s2)*(x0-x) + s2*(x1-x) ) / ( tdenom_x );
            }
            else
            {
                t2 = ( (1-s2)*(y0-y) + s2*(y1-y) ) / ( tdenom_y );
            }
            if ( !in_range( t2, 0, 1, 1e-10 ) )
                t2_valid = 0;
        }
    }

    /* Final cleanup */
    if ( t2_valid && !t_valid )
    {
        s = s2;
        t = t2;
        t_valid = t2_valid;
        t2_valid = 0;
    }

    /* Output */
    if ( t_valid )
    {
        *sout = s;
        *tout = t;
    }

    if ( t2_valid )
    {
        *s2out = s2;
        *t2out = t2;
    }

    return t_valid + t2_valid;
}

void bilerp( double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3, double s, double t, double* x, double* y )
{
    *x = t*(s*x3+(1-s)*x2) + (1-t)*(s*x1+(1-s)*x0);
    *y = t*(s*y3+(1-s)*y2) + (1-t)*(s*y1+(1-s)*y0);
}

double randrange( double range_min, double range_max )
{
    double range_width = range_max - range_min;
    double rand01 = (rand() / (double)RAND_MAX);
    return (rand01 * range_width) + range_min;
}

/* Returns number of failed trials */
int fuzzTestInvBilerp( int num_trials )
{
    int num_failed = 0;

    double x0, y0, x1, y1, x2, y2, x3, y3, x, y, s, t, s2, t2, orig_s, orig_t;
    int num_st;
    int itrial;
    for ( itrial = 0; itrial < num_trials; itrial++ )
    {
        int failed = 0;
        /* Get random positions for the corners of the quad */
        x0 = randrange( -10, 10 );
        y0 = randrange( -10, 10 );
        x1 = randrange( -10, 10 );
        y1 = randrange( -10, 10 );
        x2 = randrange( -10, 10 );
        y2 = randrange( -10, 10 );
        x3 = randrange( -10, 10 );
        y3 = randrange( -10, 10 );
        /*x0 = 0, y0 = 0, x1 = 1, y1 = 0, x2 = 0, y2 = 1, x3 = 1, y3 = 1;*/
        /* Get random s and t */
        s = randrange( 0, 1 );
        t = randrange( 0, 1 );
        orig_s = s;
        orig_t = t;
        /* bilerp to get x and y */
        bilerp( x0, y0, x1, y1, x2, y2, x3, y3, s, t, &x, &y );
        /* invert */
        num_st = inverseBilerp( x0, y0, x1, y1, x2, y2, x3, y3, x, y, &s, &t, &s2, &t2 );
        if ( num_st == 0 )
        {
            failed = 1;
        }
        else if ( num_st == 1 )
        {
            if ( !(equals( orig_s, s, 1e-5 ) && equals( orig_t, t, 1e-5 )) )
                failed = 1;
        }
        else if ( num_st == 2 )
        {
            if ( !((equals( orig_s, s , 1e-5 ) && equals( orig_t, t , 1e-5 )) ||
                   (equals( orig_s, s2, 1e-5 ) && equals( orig_t, t2, 1e-5 )) ) )
               failed = 1;
        }

        if ( failed )
        {
            num_failed++;
            printf("Failed trial %d\n", itrial);
        }
    }

    return num_failed;
}

int main( int argc, char** argv )
{
    int num_failed;
    srand( 0 );

    num_failed = fuzzTestInvBilerp( 100000000 );

    printf("%d of the tests failed\n", num_failed);
    getc(stdin);

    return 0;
}

Что ж, если p является 2D-точкой, да, вы можете легко ее получить. В этом случае S - это дробная составляющая общей ширины четырехугольника в точке T, T - также дробная составляющая общей высоты четырехугольника в точке S.

Если, однако, p является скаляром, это не обязательно возможно, потому что билинейная интерполяционная функция не обязательно монолитна.