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

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_i = (new_r * sin(theta) + center);
const double new_j = (new_r * cos(theta) + center);
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_i = center + (new_r / r) * y_affined;
const double new_j = center + (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;
const double y_translated = source_xy.y -center;
const double r = std::hypot(x_translated, y_translated);
const uint64_t new_i = static_cast<uint64_t>(center + (new_r / r) * y_translated);
const uint64_t new_j = static_cast<uint64_t>(center + (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, new_dpoint};
if (bounds.contains(new_point)) {
imageData.at<double>(i, j) = bilinearInterpolate(cpy, new_dpoint, new_dpoint);
}
}
}
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);

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

• Is there a good reason to use cv::norm(cv::Vec2d()) instead of plain old std::hypot()? May 25 '17 at 16:13
• @TobySpeight No, I just didn't know that existed as a standard function May 25 '17 at 16:18
• As it seems like you're not interested in getting a review of the context code, I marked it as a code-block-quote, indicating that it's not for review. May 25 '17 at 17:31
• @SimonForsberg That is correct, However, if I in the future would like a to come back to this same code and have a separate piece reviewed (for example, bilinear interpolation) would that be a duplicate? this is my first time posting here, so I'm not sure how to go about proper minified examples for review. I figured that I wanted to not waste to take up too much time looking at the code I didn't view as important, but after seeing what Toby posted about what he assumed I was doing, I figured it made more sense to show the rest anyway. May 25 '17 at 17:36
• As long as it is not the exact same code, it's not a duplicate. You can take a look at what you may and may not do after receiving answers.for general tips about posting follow-up questions. May 25 '17 at 19:08

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_i = center + (new_r / r) * y_affined;
const double new_j = center + (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!

• wow I did the math, I didn't think about the fact that I was just getting y_affined/r and x_affined/r with my version via cos(theta) and sin(theta). I guess I was just looking at the conversion from polar coordinates to Cartesian separately from what I was actually trying to accomplish. I think you made a mistake in your post though, I think you have the x_affined and y_affined swapped (since i corresponds to rows, and used sin in my code so op/adj or y/r) May 25 '17 at 16:31
• Yeah, it's easy to lose track of the big picture once you've broken it down into steps - not obvious how they inter-relate sometimes. Thanks for spotting the bug (that's why I said "something like" as a guard!). I've edited to fix that. May 25 '17 at 16:36
• Do comment to tell us whether this makes a measurable difference in your test (as you didn't give us enough to benchmark it ourselves). May 25 '17 at 16:38
• It does, by a lot. I'm afraid of showing all the code because the rest of the code does not appear to affect performance (well, gprof is acting weird and saying my parent function calling these takes up 96% of the run-time which I don't think is possible, I'm testing it out), and I feel that it might completely change the question, its only 106 lines of implementation though. If you think I should I'll work on transforming it to a minified example that can be independently verified. May 25 '17 at 16:41
• It's not essential - I was only justifying why I didn't present any figures myself. May 25 '17 at 16:48