Обратное билинейное отображение

При мэппинге и выполнении калибровок камер с помощью билинейного отображения иногда требуется вычислить его обратное преобразование.

Здесь приводится итеративный алгоритм вычисления обратного билинейного отображения, для случая, когда билинейное отображение bilinear(x,y) действует на единичном квадрате, и переводит его углы в точки A,B,C,D:
bilinear(x,y) = (1-x)*(1-y)*A + (1-x)y*B + x*(1-y)*C + x*y*D

//Обратное преобразование к билинейному отображению.
//Считается итеративно, методом, напоминающим деление отрезка пополам.
//P- точка, для которой надо найти обратное преобразование,
//eps - точность. Например, если единичный квадрат представляет в пикселах поле 640x640, то
//можно взять eps = 0.5 / 640 для пиксельной точности.

ofPoint bilinearInverse( ofPoint P,   float eps ) 
{
     //начальная точка
     float x = 0.5;
     float y = 0.5;
     float rad = 0.5; //радиус поиска - предполагается что destX, destY лежит в  [0,1]x[0,1]
                              //если чуть вылазит, то можно сделать 0.7, например.

    //итерационный процесс
    ofPoint p, p1;
    while ( rad > eps ) {                       //пока радиус поиска больше заданной точности
        float bestDist = 100;
        float bx = x;
        float by = y;
        for (int Y = -1; Y <= 1; Y++) {                //делаем 9 шагов сканирования, выбираем лучшее
            for (int X = -1; X <= 1; X++) {
                p.x = x + X * rad;
                p.y = y + Y * rad;
                p1 = transform( p );
                float d = ofDist( P.x, P.y, p1.x, p1.y );
                if (d < bestDist) {
                    bestDist = d;
                    bx = p.x;
                    by = p.y;
                }
            }
        }
        //перемещаем текущую точку в найденный оптимум
        x = bx;
        y = by;
        //уменьшаем радиус поиска
        rad *= 0.5;
    }
    return ofPoint( x, y );
 }