Here is another way to achieve the same effect. Though it lacks the scale part. I was trying to understand that part and implemented here as well by asking the question in the link in my comment. Anyway, here is the code:

```
static void meshgrid(InputArray _xgv, InputArray _ygv, OutputArray _X, OutputArray _Y)
{
Mat xgv = _xgv.getMat();
Mat ygv = _ygv.getMat();
cv::repeat(xgv.reshape(1,1), ygv.total(), 1, _X);
cv::repeat(ygv.reshape(1,1).t(), 1, xgv.total(), _Y);
}
// helper function (maybe that goes somehow easier)
static void meshgrid/*Test*/(const cv::Range &xgv, const cv::Range &ygv, OutputArray _X, OutputArray _Y)
{
std::vector<float> t_x, t_y;
for (int i = xgv.start; i <= xgv.end; i++) t_x.push_back(i);
for (int i = ygv.start; i <= ygv.end; i++) t_y.push_back(i);
Mat __X, __Y;
meshgrid(cv::Mat(t_x), cv::Mat(t_y), __X, __Y);
Mat(__X).copyTo(_X);
Mat(__Y).copyTo(_Y);
}
void radialDistortion(InputArray _src, OutputArray _dst, double Cx, double Cy, double coef, bool scale = true)
{
// die if distortion parameters are not correct
CV_Assert(Cx >= 0 && Cy >= 0/* && k >= 0*/);
Mat src = _src.getMat();
int ncols = src.cols;
int nrows = src.rows;
Mat xi, yi;
meshgrid(Range(1, ncols), Range(1, nrows), xi, yi);
int imid_x = Cx;
int imid_y = Cy;
Mat xt = xi - imid_x;
Mat yt = yi - imid_y;
xt = xt.t();
yt = yt.t();
if (!xt.isContinuous())
xt = xt.clone();
if (!yt.isContinuous())
yt = yt.clone();
xt = xt.reshape(0, xt.rows * xt.cols);
yt = yt.reshape(0, yt.rows * yt.cols);
Mat r, theta;
cartToPolar(xt, yt, r, theta);
Mat s = r + coef * r.mul(r.mul(r));
Mat ut, vt;
polarToCart(s, theta, ut, vt);
Mat u, v;
if (!ut.isContinuous())
ut = ut.clone();
if (!vt.isContinuous())
vt = vt.clone();
u = ut.reshape(0, ncols) + imid_x;
v = vt.reshape(0, ncols) + imid_y;
u = u.t();
v = v.t();
remap(src, _dst, u, v, CV_INTER_LINEAR, BORDER_CONSTANT);
}
int main()
{
Mat input = imread("lena.jpg"/*, CV_LOAD_IMAGE_GRAYSCALE*/); // it works with grayscale images as well ;-)
if(!input.data || input.empty())
cout << "Problem loading image!!!" << endl;
Mat output;
radialDistortion(input, output, round(static_cast<float>(input.cols) / 2), round(static_cast<float>(input.rows) / 2), 0.000003, false);
Mat canvas = Mat::zeros(input.rows, input.cols*2+10, input.type());
input.copyTo(canvas(Range::all(), Range(0, output.cols)));
output.copyTo(canvas(Range::all(), Range(output.cols+10, output.cols*2+10)));
// if it is too big to fit on the screen, then scale it down by 2, hopefully it'll fit :-)
if(canvas.cols > 1920)
{
resize(canvas, canvas, Size(canvas.cols/2, canvas.rows/2));
}
imshow("canvas", canvas);
waitKey(0);
return 0;
}
```

if you have some time you can search more about it ;-).

have a look here. I also think that I have another version of it that it can be used. I will check about it when I return home. However, there is still something that I couldn't understand and it lacks an answer, just in case you are interested to spend some time.