I have a tight loop, where I get a camera image, undistort it and also transform it according to some transformation (e.g. a perspective transform). I already figured out to use cv::remap(...) for each operation, which is already much more efficient than using plain matrix operations.

In my understanding it should be possible to combine the lookup maps into one and call remap just once in every loop iteration. Is there a canonical way to do this? I would prefer not to implement all the interpolation stuff myself.

Note: The procedure should work with differently sized maps. In my particular case the undistortion preserves the image dimensions, while the other transformation scales the image to a different size.

Code for illustration:

// 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);
}

In the case of two general mappings, there is no choice but to use the approach suggested by @MichaelBurdinov.

However, in the special case of two mappings with known inverse mappings, an alternative approach is to compute the maps manually. This manual approach is more accurate than the double remap one, since it does not involve interpolation of coordinate maps.

In practice, most of the interesting applications match this special case. It does too in your case because your first map corresponds to image undistortion (whose inverse operation is image distortion, which is associated to a well known analytical model) and your second map corresponds to a perspective transform (whose inverse can be expressed analytically).

Computing the maps manually is actually quite easy. As stated in the documentation (link) these maps contain, for each pixel in the destination image, the (x,y) coordinates where to find the appropriate intensity in the source image. The following code snippet shows how to compute the maps manually in your case:

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);

Note: H above is your perspective transform while Kshould be the camera matrix associated with the undistorted image, so it should be what in your code is called newCameraMatrix (which BTW is not an output argument of initUndistortRectifyMap). Depending on your specific data, there might also be some additional cases to handle (e.g. division by pt(2) when it might be zero, etc).

more