Opencv: исказить обратно

initUndistortRectifyMap, связанный с одним из ответов на вопрос, который вы упомянули, действительно делает то, что вы хотите. Поскольку он используется в Remap для построения полного неискаженного изображения, он дает для каждого места в целевом изображении (неискаженном) местонахождение соответствующего пикселя в искаженном изображении, чтобы они могли использовать его цвет. Так что это действительно карта f(undistorted) = distorted.

Однако использование этой карты позволит вводить только целочисленные позиции в пределах прямоугольника изображения. К счастью, документация содержит полные уравнения.

В основном это то, что у вас есть, за исключением того, что вы пропустили предварительный шаг. Вот моя версия (это С#, но должна быть такой же):

public PointF Distort(PointF point)
{
    // To relative coordinates <- this is the step you are missing.
    double x = (point.X - cx) / fx;
    double y = (point.Y - cy) / fy;

    double r2 = x*x + y*y;

    // Radial distorsion
    double xDistort = x * (1 + k1 * r2 + k2 * r2 * r2 + k3 * r2 * r2 * r2);
    double yDistort = y * (1 + k1 * r2 + k2 * r2 * r2 + k3 * r2 * r2 * r2);

    // Tangential distorsion
    xDistort = xDistort + (2 * p1 * x * y + p2 * (r2 + 2 * x * x));
    yDistort = yDistort + (p1 * (r2 + 2 * y * y) + 2 * p2 * x * y);

    // Back to absolute coordinates.
    xDistort = xDistort * fx + cx;
    yDistort = yDistort * fy + cy;

    return new PointF((float)xDistort, (float)yDistort);
}

Вы можете легко исказить свои точки, используя ProjectPoints.

cv::Mat rVec(3, 1, cv::DataType<double>::type); // Rotation vector
rVec.at<double>(0) = 0;
rVec.at<double>(1) = 0;
rVec.at<double>(2) =0;
cv::Mat tVec(3, 1, cv::DataType<double>::type); // Translation vector
tVec.at<double>(0) =0;
tVec.at<double>(1) = 0;
tVec.at<double>(2) = 0;

cv::projectPoints(points,rVec,tVec, cameraMatrix, distCoeffs,result);

PS: в opencv 3 добавили функцию для искажения.

Если вы умножите все коэффициенты искажения на -1, вы можете затем передать их в undistort или undistortPoints, и в основном вы примените обратное искажение, которое вернет искажение.

Модель камеры OCV (см. http://docs.opencv.org/2.4/modules/calib3d/doc/camera_dication_and_3d_reconstruction.html) описывает, как 3D-точка сначала сопоставляется с воображаемой идеальной координатой камеры-обскуры, а затем «искажает» координату, чтобы она моделировала изображение реальной камеры реального мира.

Используя коэффициенты искажения OpenCV (= коэффициенты искажения Брауна), можно легко рассчитать следующие 2 операции:

• Вычислить координату пикселя в исходном изображении с камеры по заданной координате пикселя в изображении без искажений (т.е. неискаженном). AFAIK для этого нет явной функции OpenCV. Но код в ответе Джоан Чарман делает именно это.
• Рассчитайте изображение без искажений на основе исходного изображения с камеры. Это можно сделать с помощью cv::undistort(....) или комбинации cv::initUndistortRectifyMap(....) и cv::remap(....).

Однако следующие 2 операции вычислительно намного сложнее:

• Вычислите координату пикселя в изображении без искажений по координате пикселя в исходном изображении с камеры. Это можно сделать с помощью cv::undistortPoints(....).
• Рассчитайте исходное изображение камеры из изображения без искажений.

Это может показаться нелогичным. Более подробное объяснение:

Для заданной координаты пикселя в изображении без искажений легко вычислить соответствующую координату в исходном изображении (т.е. «исказить» координату).

x = (u - cx) / fx; // u and v are distortion free
y = (v - cy) / fy;

rr = x*x + y*y
distortion = 1 + rr  * (k1 + rr * (k2 + rr * k3))
# I ommit the tangential parameters for clarity

u_ = fx * distortion * x + cx
v_ = fy * distortion * y + cy
// u_ and v_ are coordinates in the original camera image

Сделать это наоборот гораздо сложнее; в основном нужно было бы объединить все приведенные выше строки кода в одно большое векторное уравнение и решить его для u и v. Я думаю, что для общего случая, когда используются все 5 коэффициентов искажения, это можно сделать только численно. Что (не глядя на код), вероятно, и делает cv::undistortPoints(....).

Однако, используя коэффициенты искажения, мы можем рассчитать карту неискажения (cv::initUndistortRectifyMap(....)), которая отображает координаты изображения без искажений в исходные координаты изображения камеры. Каждая запись в карте неискажения содержит позицию пикселя (с плавающей запятой) в исходном изображении с камеры. Другими словами, карта неискажений указывает от изображения без искажений к исходному изображению с камеры. Таким образом, карта рассчитывается именно по приведенной выше формуле.

Затем карту можно применить для получения нового изображения без искажений из оригинала (cv::remap(....)). cv::undistort() делает это без явного расчета карты неискажения.

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

1. Код в вопросе почти правильный, но есть ошибка. Он использует r^4 после k3 вместо r^6. Я переписал код и успешно запустился после этого простого исправления.

xCorrected = x * (1. + k1 * r2 + k2 * r2 * r2 + k3 * r2 * r2 * r2); // Multiply r2 after k3 one more time in yCorrected too

2. Джоан Чарман сказала, что правильные вещи и уравнения позволяют искажать точку назад (не искажать, как упоминалось в комментариях под его ответом). Однако некоторые слова неверны:

// К относительным координатам ‹- это шаг, который вы пропустили

Это неправильно, так как код в этом вопросе уже использует относительные координаты! Это трюк в функциях OpenCV undistortPoints. Он имеет новую внутреннюю матрицу в качестве 6-го аргумента. Если None, то функция возвращает точки в относительных координатах. И именно поэтому исходный код, о котором идет речь, имеет этот шаг:

//Step 2 : ideal coordinates => actual coordinates
      xCorrected = xCorrected * fx + ux;
      yCorrected = yCorrected * fy + uy;

3. Также я должен сказать о неразберихе в интернет-материалах.

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

Недавно я нашел, почему. Учебник OpenCV и его документация имеют разные названия. В учебнике для уравнений используются переменные «xCorrected» и «yCorrected». Хотя в документе одни и те же вещи называются по-разному: «xDistorted» и «yDistorted».

Итак, позвольте мне разрешить путаницу: операцию искажения можно представить в виде уравнений в различных моделях искажения. Но Undistortion возможен только с помощью численного алгоритма итерации. Не существует аналитических решений для представления неискажений в виде уравнений (из-за части 6-го порядка в радиальной части и нелинейности)

Нет аналитического решения этой проблемы, если вы искажаете координаты, нет возможности вернуться назад, по крайней мере, аналитически с этой конкретной моделью. Это характерно для модели радиальной дисторсии, то, как она определена, позволяет искажать простым аналитическим способом, но не наоборот. Для этого необходимо решить многочлен 7-й степени, для которого доказано, что аналитического решения не существует.

Однако модель радиальной камеры никоим образом не является особенной или священной, это просто простое правило, которое растягивает пиксели наружу или внутрь к оптическому центру в зависимости от объектива, с которым вы делали снимок. Чем ближе к оптическому центру, тем меньше искажений получает пиксель. Существует множество других способов определения модели радиальной дисторсии, которые могут дать не только аналогичное качество дисторсии, но и обеспечить простой способ определения обратной дисторсии. Но такой путь означает, что оптимальные параметры для такой модели вам придется искать самостоятельно.

Например, в моем конкретном случае я обнаружил, что простая сигмовидная функция (смещенная и масштабированная) способна аппроксимировать существующие параметры моей радиальной модели с интегральной ошибкой MSE, меньшей или равной 1E-06, хотя сравнение между моделями кажется бессмысленным. Я не думаю, что родная радиальная модель дает лучшие значения и не должна рассматриваться как эталонная. Физическая геометрия линзы может различаться таким образом, что обе модели не могут быть представлены, и для более точного приближения геометрии линзы следует использовать подход, подобный сетке. Однако я впечатлен аппроксимированной моделью, потому что она использует только один свободный параметр и обеспечивает особенно точный результат, который заставляет меня задуматься, какая модель на самом деле лучше для работы.

Вот график исходной радиальной модели (красный) и ее сигмовидного приближения (зеленый) вверху, а также их производных (синие линии):

Итак, функция искажения/неискажения в моем случае выглядела так:

distort = (r, alpha) -> 2/(1 + exp(-alpha*r)) - 1
undistort = (d, alpha) -> -ln((d + 1)/(d - 1))/alpha

(Обратите внимание, что искажение выполняется в полярных координатах вокруг оптического центра и влияет только на расстояние от оптического центра (т.е. не на сам угол), r - расстояние от оптического центра, альфа - свободный параметр, который необходимо оценить):

Вот как выглядело искажение по сравнению с исходным радиальным искажением (зеленый цвет аппроксимирован, красный — исходное радиальное искажение)

А вот как выглядит обратное отображение пикселей, если мы возьмем обычную пиксельную сетку и попытаемся ее не исказить:

Другой способ - использовать переназначение для проецирования исправленного изображения в искаженное изображение:

img_distored = cv2.remap(img_rect, mapx, mapy, cv2.INTER_LINEAR)

mapx и mapy — это сопоставления от выпрямленных пиксельных местоположений к искаженным пиксельным местоположениям. Его можно получить, выполнив следующие шаги:

X, Y = np.meshgrid(range(w), range(h)
pnts_distorted = np.merge(X, Y).reshape(w*h, 2)
pnts_rectified = cv2.undistortPoints(pnts_distorted, cameraMatrix, distort, R=rotation, P=pose)
mapx = pnts_rectified[:,:,0]
mapy = pnts_rectified[:,:,1]

cameraMatrix, искажение, вращение, поза — это параметры, возвращаемые в функциях калибровки cv и StereoRectify.

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

В приведенном выше случае мне очень помог следующий ответ. (Руководствуясь ответом @Joan Charmant.)

distortBack(point3f undistored, point2f distored_back){
    double temp_x, temp_y;
    temp_x = undistored.x / undistored.z; 
    temp_y = undistored.y / undistored.z;
    
    double r2 = temp_x*temp_x + temp_y*temp_y;

    // Radial distorsion
    double xDistort = temp_x * (1 + dist_k1 * r2 + dist_k2 * r2 * r2 + dist_k3 * r2 * r2 * r2);
    double yDistort = temp_y * (1 + dist_k1 * r2 + dist_k2 * r2 * r2 + dist_k3 * r2 * r2 * r2);

    // Tangential distorsion
    xDistort = xDistort + (2 * dist_p1 * temp_x * temp_y + dist_p2 * (r2 + 2 * temp_x * temp_x));
    yDistort = yDistort + (dist_p1 * (r2 + 2 * temp_y * temp_y) + 2 * dist_p2 * temp_x * temp_y);

    // Back to absolute coordinates.
    distored_back.x = xDistort * camera_fx + camera_cx;
    distored_back.y = yDistort * camera_fy + camera_cy;
}

Дополнительные сведения — cv::UndistorPoints() и в opencv/doc/v4.3.0

Вот реализация, основанная на cv2.projectPoints, источник для interpolate_missing_pixels можно найти здесь.

import numpy as np
import cv2


def distort(
        image: np.ndarray,
        cam_matrix: np.ndarray,
        dist_coefs,
        **kwargs
) -> np.ndarray:
    """
    Applies lens distortion to an image.
    :param image: input image
    :param cam_matrix: camera intrinsic matrix
    :param dist_coefs must be given in opencv format
    :param kwargs passed to `interpolate_missing_pixels`.
    """
    rvec_tvec = np.zeros(3)
    h, w = image.shape[:2]
    rows, cols = np.meshgrid(np.arange(h), np.arange(w), indexing='ij')
    coords_pix_orig = np.array([cols.ravel(), rows.ravel(), np.ones(h*w)])
    coords_cam = np.linalg.inv(cam_matrix) @ coords_pix_orig
    coords_cam = coords_cam.astype(np.float32)
    coords_pix_dist, _ = cv2.projectPoints(
        coords_cam, rvec_tvec, rvec_tvec, cam_matrix, dist_coefs)
    coords_pix_dist = coords_pix_dist.astype(np.int).squeeze(1)
    in_image =\
        (coords_pix_dist[:, 0] >= 0) & (coords_pix_dist[:, 0] < w) & \
        (coords_pix_dist[:, 1] >= 0) & (coords_pix_dist[:, 1] < h)

    orig_r = coords_pix_orig[1][in_image].astype(np.int)
    orig_c = coords_pix_orig[0][in_image].astype(np.int)

    dist_r = coords_pix_dist[:, 1][in_image]
    dist_c = coords_pix_dist[:, 0][in_image]

    dist_image = np.zeros_like(image)
    dist_image[dist_r, dist_c] = image[orig_r, orig_c]

    missing_vals_mask = np.ones((h, w), dtype=np.bool)
    missing_vals_mask[dist_r, dist_c] = False

    interp_dist_image = interpolate_missing_pixels(
        dist_image, missing_vals_mask, **kwargs
    )

    return interp_dist_image