Понимание устаревшего кода: алгоритм устранения искажения радиальной линзы

Приведенный ниже метод удаления искажений объектива с камеры был написан более десяти лет назад, и я пытаюсь понять, как работает аппроксимация.

void Distortion::removeLensDistortion(const V2 &dist, V2 &ideal) const
{
 const int ITERATIONS=10;
 V2 xn;
 xn.x=(dist.x - lensParam.centerX)/lensParam.fX;
 xn.y=(dist.y - lensParam.centerY)/lensParam.fY;

 V2 x=xn;
 for (int i=0;i<ITERATIONS;i++) {
   double r2 = Utilities::square(x.x)+Utilities::square(x.y);
   double r4 = Utilities::square(r2);
   double rad=1+lensParam.kc1 * r2 + lensParam.kc2 * r4;
   x.x/=rad;
   x.y/=rad;
 }
 ideal.x=x.x*lensParam.fX+lensParam.centerX;
 ideal.y=x.y*lensParam.fY+lensParam.centerY;
}

Как напоминание:

• линзаParam.centerX и линзаParam.centerY является главной точкой
• LenParam.fX и LenParam.fY — фокусное расстояние в пикселях.
• линзыParam.kc1 и линзыParam.kc2 — первые два коэффициента радиальной дисторсии. Это k_1 и k_2 в приведенной ниже формуле.

Формула для добавления дисторсии объектива с учетом первых двух параметров радиальной дисторсии выглядит следующим образом:

x_distorted = x_undistorted * (1+k_1 * r² + k_2 * r^4) 
y_distorted = y_undistorted * (1+k_1 * r² + k_2 * r^4)

где r²=(x_undistorted)²+(y_undistorted)² и r^4=(r²)²

В приведенном выше коде термин (1+k_1 * r² + k_2 * r^4) вычисляется и сохраняется в переменной rad, а искаженный x делится на rad в каждой из десяти итераций. Все камеры, которые мы используем, имеют подушкообразную дисторсию (так что k_1<0)

Вопрос в том, как этот алгоритм аппроксимирует неискаженные точки изображения? Знаете ли вы, есть ли какая-либо статья, в которой предлагается этот алгоритм?

Неискажение opencv может быть немного похоже, поэтому эта ссылка может быть полезной, но это не совсем хотя то же самое.