Backward transformation implementation for a barrel eye lens aberration in c++ 14 (using opencv)

Question

I've made the following backwards transformation from a "transformed" source x,y point, to the resulting index of the actual image that I have, so I can avoid black spots appearing in my barrel eye lens image.

I've used the formula focal_length * arctan(radius, focal_length) for image undistortion. Here is the code I've used (its part of another class, but the rest of the class is not important):

double BarrelAbberationClass::toDistortedRadius(const double r) const {
    return m_focal_length * atan2(r, m_focal_length);
}

cv::Vec2d BarrelAbberationClass::transformBackward(const cv::Point &source_xy,
                                               const cv::Mat &affine, const cv::Vec2d &center) {
    const double x_affined = source_xy.x * affine.at<double>(0, 0) + affine.at<double>(0, 2);
    const double y_affined = source_xy.y * affine.at<double>(1, 1) + affine.at<double>(1, 2);
    const double r = cv::norm(cv::Vec2d(x_affined, y_affined));
    const double theta = atan2(y_affined, x_affined);
    const double new_r = toUnDistortedRadius(r);
    const double new_i = (new_r * sin(theta) + center[1]);
    const double new_j = (new_r * cos(theta) + center[0]);
    return cv::Vec2d(new_j, new_i);
}

Basically, given a source x,y point, affine matrix for the transformation account for the center of my image actually being width/2 and height/2, and the visual bounds I would even be able to grab from in the non transformed image, then I warp the x and y affined points to the correct undistorted plane and return the calculated row and column i, j that would correspond to points in the original image (these must be floating point for use later in bilinear interpolation of positions in between others)

I should also note that the code with affine.at<double>(...) were originally matrix operations. I saw in the profiler that this was slowing doing my code. Removing the matrix creation routines caused by matrix operations here increased my speed by 4.

When profiling this code (I'm forced to use Windows + MingW for this, so I don't have many convenient profiling options; I've been using gprof) it looks like tan, cos, sin, and atan2 take up the vast majority of the time at -O3, with atan2 cos and sin each taking about 20% of the time.

Here is the full implementation to give context:

double BarrelAbberationClass::toDistortedRadius(const double r) const {
    return m_focal_length * atan2(r, m_focal_length);
}

cv::Vec2d BarrelAbberationClass::transformBackward(const cv::Point &source_xy,
                                               const cv::Mat &affine, const cv::Vec2d &center) {
    const double x_affined = source_xy.x * affine.at<double>(0, 0) + affine.at<double>(0, 2);
    const double y_affined = source_xy.y * affine.at<double>(1, 1) + affine.at<double>(1, 2);
    const double r = std::hypot(x_affined, y_affined);
    const double new_r = toUnDistortedRadius(r);
    const double new_i = center[1] + (new_r / r) * y_affined;
    const double new_j = center[0] + (new_r / r) * x_affined;
    return cv::Vec2d(new_j, new_i);
}


double bilinearInterpolate(const cv::Mat &cv_source, const float x, const float y) {
    const int px = static_cast<int>(x);
    const int py = static_cast<int>(y);

    const double p1 = cv_source.at<double>(py, px);
    const double p2 = cv_source.at<double>(py, px + 1);
    const double p3 = cv_source.at<double>(py + 1, px);
    const double p4 = cv_source.at<double>(py + 1, px + 1);

    const float fx = x - px;
    const float fy = y - py;
    const float fx1 = 1.0f - fx;
    const float fy1 = 1.0f - fy;

    const float w1 = fx1 * fy1;
    const float w2 = fx * fy1;
    const float w3 = fx1 * fy;
    const float w4 = fx * fy;

    return p1 * w1 + p2 * w2 + p3 * w3 + p4 * w4;
}

double BarrelAbberationClass::toUnDistortedRadius(const double r) const {
    return m_focal_length * tan(r / m_focal_length);
}

cv::Point BarrelAbberationClass::transformForward(const cv::Point &source_xy,
                                              const cv::Vec2d &center) {
    const double x_translated = source_xy.x -center[0];
    const double y_translated = source_xy.y -center[1];
    const double r = std::hypot(x_translated, y_translated);
    const double new_r = toDistortedRadius(r);
    const uint64_t new_i = static_cast<uint64_t>(center[1] + (new_r / r) * y_translated);
    const uint64_t new_j = static_cast<uint64_t>(center[0] + (new_r / r) * x_translated);
    return cv::Point(new_j, new_i);
}

cv::Mat BarrelAbberationClass::getScaleMatrix(const cv::Rect &bounds,
                                          const cv::Vec2d &center) {
    const cv::Point new_tl = transformForward(bounds.tl(), center);
    const cv::Point new_br = transformForward(bounds.br(), center);
    const cv::Point diff = new_br - new_tl;
    const double scale_x = diff.x / static_cast<double>(bounds.width);
    const double scale_y = diff.y / static_cast<double>(bounds.height);
    const double scale = scale_x > scale_y ? scale_x : scale_y;
    cv::Mat scale_matrix = (cv::Mat_<double>(3, 3) << scale, 0, 0,
            0, scale, 0,
            0, 0, 1);
    return scale_matrix;
}

cv::Mat& BarrelAbberationClass::calculateAberation(cv::Mat& imageData) {
    const double center_x = imageData.size().width / 2;
    const double center_y = imageData.size().height / 2;
    const cv::Vec2d center(center_x, center_y);
    const cv::Mat translate = (cv::Mat_<double>(3, 3) << 1, 0, -center_x,
                                                         0, 1, -center_y,
                                                         0, 0, 1);
    const cv::Mat cpy = imageData.clone();
    imageData.setTo(cv::Scalar(0));
    const cv::Rect bounds(cv::Point(), cpy.size());
    const cv::Mat scale = getScaleMatrix(bounds, center);
    const cv::Mat affine = scale * translate;

    for (uint64_t i = 0; i < cpy.rows; ++i) {
        for (uint64_t j = 0; j < cpy.cols; ++j) {
            const cv::Vec2d new_dpoint = transformBackward(cv::Point(j, i), affine, center);
            const cv::Point new_point = {new_dpoint[0], new_dpoint[1]};
            if (bounds.contains(new_point)) {
                imageData.at<double>(i, j) = bilinearInterpolate(cpy, new_dpoint[0], new_dpoint[1]);
            }
        }
    }
    return imageData;
}

Here is the .h file

class BarrelAbberationClass{

protected:
    const double m_focal_length;

public:

    BarrelDegredation(const double focal_length) : m_focal_length(focal_length) {};

    cv::Mat& calculateAberation(cv::Mat& imageData);

    double toDistortedRadius(const double r) const;

    double toUnDistortedRadius(const double r) const;

    cv::Point transformForward(const cv::Point &source_xy, const cv::Vec2d &center);

    cv::Mat getScaleMatrix(const cv::Rect &bounds, const cv::Vec2d &center);

    cv::Vec2d transformBackward(const cv::Point &source_xy,
                                const cv::Mat &affine, const cv::Vec2d &center);
}

Answer 1

It might make no difference to performance, but it may be a little better to use std::hypot() rather than create a temporary cv::Vec2d.

You can avoid working with the angle theta, simply by using the ratio new_r/r to scale the distance from centre. Something like:

const double r = std::hypot(x_affined, y_affined);
const double new_r = toUnDistortedRadius(r);
const double new_i = center[1] + (new_r / r) * y_affined;
const double new_j = center[0] + (new_r / r) * x_affined;

I've parenthesized new_r / r just in case your compiler needs extra help to hoist this common expression.

---

I'm guessing that you're passing this function to a general-purpose transform function in OpenCV. If you have control over the order in which the output pixels are iterated, and if the affine transformation has no shear components (just translation and rotation), you might be able to reduce the calculation of r and new_r almost eightfold by observing that (relative to the centre position), r is the same for ±x,±y and ±y,±x. This will hit your locality of reference, though, so may not be as significant as it sounds!