如何将两个remap()操作合并为一个？

Question

我有一个紧凑的循环，其中我获取摄像头图像，对其进行去畸变和一些变换（例如透视变换），并利用cv::remap(...)完成每个操作，这比使用纯矩阵运算要高效得多。

据我了解，可以将查找映射组合成一个，并在每个循环迭代中仅调用一次remap。是否有一种通用的方法来实现此目的？我不想自己实现所有的插值处理。

注意：该过程应适用于不同大小的映射。在我的特定情况下，去畸变保留图像尺寸，而其他变换则将图像缩放为不同的大小。

示例代码：

// input arguments
const cv::Mat_<math::flt> intrinsic  = getIntrinsic();
const cv::Mat_<math::flt> distortion = getDistortion();
const cv::Mat mNewCameraMatrix = cv::getOptimalNewCameraMatrix(intrinsic, distortion, myImageSize, 0);

// output arguments
cv::Mat undistortMapX;
cv::Mat undistortMapY;

// computes undistortion maps
cv::initUndistortRectifyMap(intrinsic, distortion, cv::Mat(),
                            newCameraMatrix, myImageSize, CV_16SC2,
                            undistortMapX, undistortMapY);

// computes undistortion maps
// ...computation of mapX and mapY omitted
cv::convertMaps(mapX, mapY, skewMapX, skewMapY, CV_16SC2);

for(;;) {
    cv::Mat originalImage = getNewImage();

    cv::Mat undistortedImage;
    cv::remap(originalImage, undistortedImage, undistortMapX, undistortMapY, cv::INTER_LINEAR);

    cv::Mat skewedImage;
    cv::remap(undistortedImage, skewedImage, skewMapX, skewMapY, cv::INTER_LINEAR);

    outputImage(skewedImage);
}

- Dimitri Schachmann

Answer 1

您可以在undistortMapX和undistortMapY上应用remap。

cv::remap(undistortMapX, undistrtSkewX, skewMapX, skewMapY, cv::INTER_LINEAR);
cv::remap(undistortMapY, undistrtSkewY, skewMapX, skewMapY, cv::INTER_LINEAR);

那么你可以使用：

cv::remap(originalImage , skewedImage, undistrtSkewX, undistrtSkewY, cv::INTER_LINEAR);

它能正常运行的原因是skewMaps和undistortMaps是图像中坐标的数组，因此它应该类似于定位位置的定位...

编辑（回复评论）：

我想我需要做一些澄清。remap()函数从旧图像的像素计算新图像中的像素。在线性插值的情况下，新图像中的每个像素是来自旧图像的4个像素的加权平均值。权重根据提供的映射值从像素到像素有所不同。如果该值接近整数，则大部分权重将来自单个像素。结果新图像将与原始图像一样锐利。另一方面，如果该值远离整数（即整数+0.5），则权重相似。这将产生平滑效果。要了解我所说的内容，请查看无畸变图像。您将看到图像的某些部分比其他部分更锐利/平滑。

现在回到关于将两个remap操作合并为一个时发生了什么的解释。组合映射中的坐标是正确的，即从原始图像的正确4个像素按正确的权重计算skewedImage的像素。但它与两次remap操作的结果不完全相同。无畸变图像中的每个像素都是来自原始图像的4个像素的加权平均值。这意味着skewedImage的每个像素将是来自orginalImage的9-16个像素的加权平均值。结论：使用单个remap()不可能给出与两个remap()使用完全相同的结果。

关于哪种可能的图像（单一的remap() vs 双重的remap()）更好的讨论是相当复杂的。通常最好尽可能少地进行插值，因为每次插值会引入不同的伪影。特别是如果伪影在图像中不均匀（某些区域比其他区域更光滑）。在某些情况下，这些伪影可能对图像产生良好的视觉效果-如减少一些抖动。但如果这正是您想要的，那么您可以以更便宜和更一致的方式实现此目标。例如，在重新映射之前平滑原始图像。

- Michael Burdinov

---

谢谢！我测试了一下，总体的去畸变和图像倾斜看起来不错，但插值似乎没有正常工作。之前我有平滑的对角线，现在它们在小范围内有点锯齿状。 - Dimitri Schachmann

不客气。你应用的第一个remap（去畸变）平滑了原始图像的一部分。在矫正畸变的情况下，这通常看起来像是平滑和未平滑区域的条纹。如果第二个remap（倾斜）从平滑区域取像素，则应该会产生相对平滑的对角线。如果它从不太平滑的区域取像素，则生成的线条会有些锯齿状。如果这些锯齿状的线条是一个问题，那么无论如何，在重新映射之前都应该平滑originalImage。 - Michael Burdinov

实际上，我不确定您所说的“通过第一次重新映射平滑”的意思是什么。我看到的是，生成的图像与我的两个重新映射的原始方法不同。在大尺度上相同，但是当您仔细观察时会有抖动。如果您希望，我可以稍后发布一些屏幕截图。我绝对需要的是，结果与在每个循环迭代中使用两个“重新映射”完全相同。 - Dimitri Schachmann

很棒的答案——效果不错。越界像素是黑色和其他颜色的混合。也许需要特别处理越界像素？ - wcochran

Answer 2

对于两个一般映射，只能使用@MichaelBurdinov建议的方法。

然而，在已知逆映射的情况下，两个映射的特殊情况可以采用手动计算地图的替代方法。这种手动方法比双重重新映射方法更精确，因为它不涉及坐标映射的插值。

在实践中，大多数有趣的应用程序都符合这种特殊情况。在您的情况下，第一个映射对应于图像失真校正（其逆操作是图像畸变，与众所周知的分析模型相关），而第二个映射对应于透视变换（其逆可以用解析法表达）。

手动计算地图实际上很容易。如文档所述（链接），这些地图包含目标图像中每个像素的(x,y)坐标，用于在源图像中找到相应的强度。以下代码片段展示了如何在您的情况下手动计算地图：

int dst_width=...,dst_height=...;           // Initialize the size of the output image
cv::Mat Hinv=H.inv(), Kinv=K.inv();         // Precompute the inverse perspective matrix and the inverse camera matrix
cv::Mat map_undist_warped_x32f(dst_height,dst_width,CV_32F);    // Allocate the x map to the correct size (n.b. the data type used is float)
cv::Mat map_undist_warped_y32f(dst_height,dst_width,CV_32F);    // Allocate the y map to the correct size (n.b. the data type used is float)
// Loop on the rows of the output image
for(int y=0; y<dst_height; ++y) {
    std::vector<cv::Point3f> pts_undist_norm(dst_width);
    // For each pixel on the current row, first use the inverse perspective mapping, then multiply by the
    // inverse camera matrix (i.e. map from pixels to normalized coordinates to prepare use of projectPoints function)
    for(int x=0; x<dst_width; ++x) {
        cv::Mat_<float> pt(3,1); pt << x,y,1;
        pt = Kinv*Hinv*pt;
        pts_undist_norm[x].x = pt(0)/pt(2);
        pts_undist_norm[x].y = pt(1)/pt(2);
        pts_undist_norm[x].z = 1;
    }
    // For each pixel on the current row, compose with the inverse undistortion mapping (i.e. the distortion
    // mapping) using projectPoints function
    std::vector<cv::Point2f> pts_dist;
    cv::projectPoints(pts_undist_norm,cv::Mat::zeros(3,1,CV_32F),cv::Mat::zeros(3,1,CV_32F),intrinsic,distortion,pts_dist);
    // Store the result in the appropriate pixel of the output maps
    for(int x=0; x<dst_width; ++x) {
        map_undist_warped_x32f.at<float>(y,x) = pts_dist[x].x;
        map_undist_warped_y32f.at<float>(y,x) = pts_dist[x].y;
    }
}
// Finally, convert the float maps to signed-integer maps for best efficiency of the remap function
cv::Mat map_undist_warped_x16s,map_undist_warped_y16s;
cv::convertMaps(map_undist_warped_x32f,map_undist_warped_y32f,map_undist_warped_x16s,map_undist_warped_y16s,CV_16SC2);

注意：H是您的透视变换，而K应该是与无畸变图像相关联的相机矩阵，因此它应该是您的代码中称为newCameraMatrix的内容（顺便说一句，这不是initUndistortRectifyMap的输出参数）。根据您特定的数据情况，可能还有其他一些需要处理的情况（例如，当pt(2)可能为零时除法等）。

- BConic

---

很棒的答案，算法对我来说也是如期工作的。img2.cols 应该改为 dst_width。我还将更正我的代码片段以反映 newCameraMatrix 的来源。 - Dimitri Schachmann

我花了一些时间试图找出为什么我的结果与Michael Burdinov的答案不同，然后我意识到我的地图是由鱼眼库（cv::fisheye::initUndistortRectifyMap()）生成的，所以如果您也使用鱼眼，则需要使用来自鱼眼库的projectPoints函数才能使此答案适用于您：cv::fisheye::projectPoints()。 - BourbonCreams

我个人需要将点的声明从<float>更改为<double>。 - BourbonCreams

Answer 3

当我在寻找如何在Python中结合去畸变和投影转换时，我发现了这个问题，但是没有直接的Python答案。

这里是BConic答案的Python直接转换。

import numpy as np
import cv2

dst_width = ...
dst_height = ...

h_inv = np.linalg.inv(h)
k_inv = np.linalg.inv(new_camera_matrix)

map_x = np.zeros((dst_height, dst_width), dtype=np.float32)
map_y = np.zeros((dst_height, dst_width), dtype=np.float32)

for y in range(dst_height):
    pts_undist_norm = np.zeros((dst_width, 3, 1))
    for x in range(dst_width):
        pt = np.array([x, y, 1]).reshape(3,1)
        pt2 = k_inv @ h_inv @ pt
        pts_undist_norm[x][0] = pt2[0]/pt2[2] 
        pts_undist_norm[x][1] = pt2[1]/pt2[2]
        pts_undist_norm[x][2] = 1

    r_vec = np.zeros((3,1))
    t_vec = np.zeros((3,1))
    pts_dist, _ = cv2.projectPoints(pts_undist_norm, r_vec, t_vec, intrinsic, distortion)

    pts_dist = pts_dist.squeeze()
    for x2 in range(dst_width):
        map_x[y][x2] = pts_dist[x2][0]
        map_y[y][x2] = pts_dist[x2][1]

# using CV_16SC2 introduced substantial image artifacts for me
map_x_final, map_y_final = cv2.convertMaps(map_x, map_y, cv2.CV_32FC1, cv2.CV_32FC1) 

显然这种方法非常缓慢，因为它使用双重循环并遍历每个像素，所以您可以使用numpy更快地完成此操作。您应该能够在C++中执行类似的操作，消除for循环并进行单矩阵乘法。

import numpy as np
import cv2

dst_width = ...
dst_height = ...

h_inv = np.linalg.inv(h)
k_inv = np.linalg.inv(new_camera_matrix)

m_grid = np.mgrid[0:dst_width, 0:dst_height].reshape(2, dst_height*dst_width)
m_grid = np.insert(m_grid, 2, 1, axis=0)

m_grid_result = k_inv @ h_inv @ m_grid
pts_undist_norm = m_grid_result[:2, :] / m_grid_result[2, :]
pts_undist_norm = np.insert(pts_undist_norm, 2, 1, axis=0)

r_vec = np.zeros((3,1))
t_vec = np.zeros((3,1))
pts_dist = cv2.projectPoints(pts_undist_norm, r_vec, t_vec, intrinsic, distortion)
pts_dist = pts_dist.squeeze().astype(np.float32)

map_x = pts_dist[:, 0].reshape(dst_width, dst_height).swapaxes(0,1)
map_y = pts_dist[:, 1].reshape(dst_width, dst_height).swapaxes(0,1)

# using CV_16SC2 introduced substantial image artifacts for me
map_x_final, map_y_final = cv2.convertMaps(map_x, map_y, cv2.CV_32FC1, cv2.CV_32FC1) 

这个numpy实现大约比第一种方法快25-75倍。

- Jarod Hanko

Answer 4

我遇到了同样的问题。我尝试实现AldurDisciple的答案。不是在循环中计算变换，而是使用一个带有mat.at <Vec2f>(x,y)=Vec2f(x,y)的矩阵，并将其应用于perspectiveTransform。将“1”的第三个通道添加到结果矩阵中，并应用projectPoints。 这是我的代码：

Mat xy(2000, 2500, CV_32FC2);
float *pxy = (float*)xy.data;
for (int y = 0; y < 2000; y++)
    for (int x = 0; x < 2500; x++)
    {
        *pxy++ = x;
        *pxy++ = y;
    }

// perspective transformation of coordinates of destination image,
// which generates the map from destination image to norm points
Mat pts_undist_norm(2000, 2500, CV_32FC2);
Mat matPerspective =transRot3x3;
perspectiveTransform(xy, pts_undist_norm, matPerspective);

//add 3rd channel of 1
vector<Mat> channels;
split(pts_undist_norm, channels);
Mat channel3(2000, 2500, CV_32FC1, cv::Scalar(float(1.0)));
channels.push_back(channel3);
Mat pts_undist_norm_3D(2000, 2500, CV_32FC3);
merge(channels, pts_undist_norm_3D);

//projectPoints to extend the map from norm points back to the original captured image  
pts_undist_norm_3D = pts_undist_norm_3D.reshape(0, 5000000);
Mat pts_dist(5000000, 1, CV_32FC2);
projectPoints(pts_undist_norm_3D, Mat::zeros(3, 1, CV_64F), Mat::zeros(3, 1, CV_64F), intrinsic, distCoeffs, pts_dist);
Mat maps[2];
pts_dist = pts_dist.reshape(0, 2000);
split(pts_dist, maps);

// apply map
remap(originalImage, skewedImage, maps[0], maps[1], INTER_LINEAR);

用于映射到规范点的变换矩阵与AldurDisciple答案中使用的略有不同。transRot3x3由calibrateCamera生成的tvec和rvec组成。

double transData[] = { 0, 0, tvecs[0].at<double>(0), 0, 0, 
tvecs[0].at<double>(1), 0, 0,  tvecs[0].at<double>(2) };
Mat translate3x3(3, 3, CV_64F, transData);
Mat rotation3x3;
Rodrigues(rvecs[0], rotation3x3);

Mat transRot3x3(3, 3, CV_64F);
rotation3x3.col(0).copyTo(transRot3x3.col(0));
rotation3x3.col(1).copyTo(transRot3x3.col(1));
translate3x3.col(2).copyTo(transRot3x3.col(2));

新增：

我意识到，如果只需要最终的地图，为什么不直接使用projectPoints将其投影到一个带有mat.at(x,y)=Vec2f(x,y,0)的矩阵中。

//generate a 3-channel mat with each entry containing it's own coordinates
Mat xyz(2000, 2500, CV_32FC3);
float *pxyz = (float*)xyz.data;
for (int y = 0; y < 2000; y++)
    for (int x = 0; x < 2500; x++)
    {
        *pxyz++ = x;
        *pxyz++ = y;
        *pxyz++ = 0;
    }

// project coordinates of destination image,
// which generates the map from destination image to source image directly
xyz=xyz.reshape(0, 5000000);
Mat pts_dist(5000000, 1, CV_32FC2);
projectPoints(xyz, rvecs[0], tvecs[0], intrinsic, distCoeffs, pts_dist);
Mat maps[2];
pts_dist = pts_dist.reshape(0, 2000);
split(pts_dist, maps);

//apply map
remap(originalImage, skewedImage, maps[0], maps[1], INTER_LINEAR);

- ricecakebear