这篇文章主要讲解桶形畸变矫正，首先介绍一下相关背景。

Radial Distortion

尽管畸变可能是没有规律或遵循某种模式的，但一般情况下我们遇到的最多的情况是畸变呈放射状且是对称的，其主要产生于相机镜头的畸变。这种呈放射状的、对称的畸变，即radial distortion主要可以分为两种，barrel distortion（桶形畸变）以及pincushion distortion（枕型畸变）。

它们的表现如下图所示，

畸变种类	效果	说明
桶形畸变		桶形畸变中，从图像的中心点到图像的边缘点图像的相对大小逐渐变小，看起来就像把一张正常大小的图片包裹在一个球上一样。鱼眼相机，就是利用这种类型的畸变将一个无限宽的物体平面映射成有限的图像区域来呈现出一种半球形状的效果的。当使用变焦镜头时如果把焦距调至镜头的中段就会出现桶形畸变，此外使用广角镜头时桶形畸变的效果最明显。
枕型畸变		枕型畸变中，从图像的中心点到图像的边缘点图像的相对大小逐渐变大，看起来的效果就像是把一张正常大小的图片贴在球的内壁，像一个枕头。
胡子型畸变		胡子型畸变是以上两种畸变的混合，有时也叫做混合型畸变。从中心向边缘扩散的过程中，畸变种类也从桶形畸变变成了枕型畸变。正因为在图片的上半部分中的直线变得像胡子形状一样，所以我们称其为“胡子型畸变”(Mustache distortion)。

从数学上来说，桶形畸变和枕型畸变都是二次的，意味着畸变程度与点到中心点之间的距离的二次方成正比。在胡子型畸变中，四次方项是显著的，是因为在中心处，桶形畸变（二次方）占主要影响，而在边缘处四次方的枕型畸变占上方。

Software Correction

低阶放射状畸变可以用Brown的矫正模型（https://web.archive.org/web/20180312205006/https://www.asprs.org/wp-content/uploads/pers/1966journal/may/1966_may_444-462.pdf）进行矫正，这个模型也叫做Brown-Conrady模型（基于Conrady早期的一些工作，http://adsabs.harvard.edu/full/1919MNRAS..79..384C）。这个Brown-Conrady模型既可以用于桶形畸变的矫正，也可以用于切线畸变的矫正（偏心矫正，discentering distortion）。

$x_u = x_d + (x_d - x_c)(K_1r^2 + K_2r^4 + ...) + (P_1(r^2 + 2(x_d - x_c)^2) + 2P_2(x_d - x_c)(y_d - y_c)) (1 + P_3r^2 + P_4r^4 + ...)$

$y_u = y_d + (y_d - y_c)(K_1r^2 + K_2r^4 + ...) + (P_2(r^2 + 2(y_d - y_c)^2) + 2P_1(x_d - x_c)(y_d - y_c)) (1 + P_3r^2 + P_4r^4 + ...)$

参数解释，

参数	含义
$(x_d, y_d)$	畸变图上的点坐标，d代表distorted
$(x_u, y_u)$	矫正图上的点坐标，u代表undistored
$(x_c, y_c)$	畸变图的中心点坐标
$K_n$	$n^{th}$ 放射畸变系数
$P_n$	$n^{th}$ 切线畸变系数
$r$	$\sqrt{(x_d - x_c)^2 + (y_d - y_c)^2}$

其中，桶形畸变的参数 $K_1$ 是负的，而枕型畸变的是正数。而胡子型畸变的放射几何序列是非单调的，即对于一些 $r$ 而言，序列的符号是变化的。在给放射畸变建模的过程中，除法模型（division model，https://ieeexplore.ieee.org/document/990465）比Brown-Conrady的偶数阶多项式模型（http://www.cs.ait.ac.th/~mdailey/papers/Bukhari-RadialDistortion.pdf）更准确。

$x_u = x_c + \frac{x_d - x_c}{1 + K_1r^2 + K_2r^4}$

$y_u = y_c + \frac{y_d - y_c}{1 + K_1r^2 + K_2r^4}$

此处的参数定义同上面那个公式相同。对于放射状畸变而言，除法模型比Brown-Conrady模型更受欢迎，主要是因为它可以使用更少的参数达到更好的效果。使用除法模型，对大多数相机而言一个调整项就已经足够了。

Python代码实现

from PIL import Image
import numpy as np
import imageio
import math
import cv2
import timedef read_image_return_nparray(fname):"""open image use Image and convert it into np.array:param fname: image file name:return: np.array(image)"""image = Image.open(fname)image_array = np.array(image)return image_arraydef read_color_image_as_gray_return_nparray(fname):image = cv2.imread(fname)original_gray_image = cv2.cvtColor(image, cv2.COLOR_RGB2GRAY)return original_gray_imagedef get_image_height_and_width_and_channels(image_array):return image_array.shape[1], image_array.shape[0], image_array.shape[2]def barrel_correction(distorted_image, image_width, image_height, paramA, paramB, paramC, paramD):distorted_image_copy = distorted_imaged = int(min(image_width, image_height) / 2)center_x = (image_width - 1) / 2.0center_y = (image_height - 1) / 2.0for x in range(image_width):for y in range(image_height):delta_x = (x - center_x) / ddelta_y = (y - center_y) / ddst_r = np.sqrt(delta_x * delta_x + delta_y * delta_y)src_r = (paramA * dst_r * dst_r * dst_r + paramB * dst_r * dst_r + paramC * dst_r + paramD) * dst_rfactor = dst_r / src_rsrc_xd = center_x + (delta_x * factor * d)src_yd = center_y + (delta_y * factor * d)src_x = int(src_xd)src_y = int(src_yd)if 0 <= src_x < image_width and 0 <= src_y < image_height:distorted_image[y][x] = distorted_image_copy[src_y][src_x]return distorted_imagedef barrel_correction2(src_image, param_a, param_b, param_c, param_d):xDim = get_image_height_and_width_and_channels(src_image)[0]yDim = get_image_height_and_width_and_channels(src_image)[1]zDim = get_image_height_and_width_and_channels(src_image)[2]dest_image = np.zeros_like(src_image)xcen = (xDim - 1.0) / 2.0ycen = (yDim - 1.0) / 2.0normDist = min(xcen, ycen)imageMin = np.min(src_image)dest_image.fill(imageMin)for k in range(zDim):for j in range(yDim):yoff = (j - ycen) / normDistfor i in range(xDim):xoff = (i - xcen) / normDistrdest2 = xoff * xoff + yoff * yoffrdest = math.sqrt(rdest2)rdest3 = rdest2 * rdestrdest4 = rdest2 * rdest2rsrc = param_a * rdest4 + param_b * rdest3 + param_c * rdest2 + param_d * rdestrsrc = normDist * rsrcang = math.atan2(yoff, xoff)xSrc = xcen + (rsrc * math.cos(ang))ySrc = ycen + (rsrc * math.sin(ang))if 0 <= xSrc < xDim - 1 and 0 <= ySrc < yDim - 1:xBase = int(math.floor(xSrc))delX = float(xSrc - xBase)yBase = int(math.floor(ySrc))delY = float(ySrc - yBase)dest_image[j][i][k] = int((1 - delX) * (1 - delY) * src_image[yBase][xBase][k])if xSrc < (xDim - 1):dest_image[j][i][k] += int(delX * (1 - delY) * src_image[yBase][xBase + 1][k])if ySrc < (yDim - 1):dest_image[j][i][k] += int((1 - delX) * delY * src_image[yBase + 1][xBase][k])if (xSrc < (xDim - 1)) and (ySrc < (yDim - 1)):dest_image[j][i][k] += int(delX * delY * src_image[yBase + 1][xBase + 1][k])return dest_imagedef save_image_array(image_array, fname):imageio.imwrite(fname, image_array)if __name__ == "__main__":start = time.time()print("start!!! " + time.ctime())original_image = read_image_return_nparray("/Users/pengyuyan/Desktop/singleFrame/images/singleImages/3.png")# original_grey_image = read_color_image_as_gray_return_nparray(#     "/Users/pengyuyan/Desktop/singleFrame/images/singleImages/4.png")h, w, z = get_image_height_and_width_and_channels(original_image)paramA = 0.08  # affects only the outermost pixels of the imageparamB = -0.36  # most cases only require b optimizationparamC = 0.0  # most uniform correctionparamD = 1.0 - paramA - paramB - paramC  # describes the linear scaling of the imagecorrected_original_image = barrel_correction2(original_image, paramA, paramB, paramC, paramD)# colored_result = cv2.cvtColor(corrected_original_image, cv2.COLOR_GRAY2RGB)save_image_array(corrected_original_image,"/Users/pengyuyan/Desktop/singleFrame/images/mytest_result_images/barrel3_result_python.png")end = time.time()print("end!!! " + time.ctime())print("duration: {}".format(end - start))

其中，paramA、paramB、paramC以及paramD需要根据自己的实际图片进行调整。这里给出一张实验图片及其对应的结果。

提升代码运行效率

以上代码的运行效率是很慢的，在经过同学的帮助之后，先是引入了numpy的一些向量运算，把速度提升到了1.5s左右，然后因为我们的眼镜参数是固定的，因此可以考虑把convert map保存成文件，然后只要通过文件读取map之后用cv2.remap函数，这样把速度提升到了每张图片100ms左右。（convert map的具体思想来自于stackoverflow，https://stackoverflow.com/questions/46520123/how-do-i-use-opencvs-remap-function）

以下是我的代码。首先第一个是获取图片之间对应关系的代码。

import numpy as np
import cv2def read_image(fname):image = cv2.imread(fname)return imagedef get_image_height_and_width(image):return image.shape[:2]def get_mapping_relation(source_image, a, b, c, d, mapping_fname):height, width = get_image_height_and_width(source_image)x, y = np.meshgrid(range(width), range(height))x = x.reshape(-1)y = y.reshape(-1)location_of_source_image = np.stack([x, y], 1)center_x = x.mean()center_y = y.mean()center = np.array([center_x, center_y])norm = np.mean(center)dist = np.sqrt(((location_of_source_image - center) ** 2).sum(1))r = np.sqrt(((x - center_x) / norm) ** 2 + ((y - center_y) / norm) ** 2)rdest = (a * r ** 4 + b * r ** 3 + c * r ** 2 + d * r) * normtarget_x = rdest / dist * (x - center_x) + center_xtarget_y = rdest / dist * (y - center_y) + center_ylocation_of_dest_image = np.stack([target_x, target_y], 1)save_array_as_npy(mapping_fname, location_of_dest_image)def save_array_as_npy(fname, array):np.save(fname, array)if __name__ == "__main__":# 1. set some necessary paramsparam_a = 0.08  # affects only the outermost pixels of the imageparam_b = -0.37  # most cases only require b optimizationparam_c = 0.0  # most uniform correctionparam_d = 1.0 - param_a - param_b - param_c  # describes the linear scaling of the imagemapping_where_to_save = "/Users/pengyuyan/Desktop/singleFrame/images/params/location_of_dest_image.npy"# 2. read a imagecolor = read_image("/Users/pengyuyan/Desktop/singleFrame/images/singleImages/4.png")# 3. get mapping and save it as npyget_mapping_relation(color, param_a, param_b, param_c, param_d, mapping_where_to_save)

第二个是把对应关系直接取出来用的代码。

import numpy as np
import cv2
import timeif __name__ == "__main__":start = time.time()print("start!!! " + time.ctime())dest_array = np.load("/Users/pengyuyan/Desktop/singleFrame/images/params/location_of_dest_image.npy")dest_array = dest_array.reshape((-1, 2))for i in range(1, 5):original_image_name = "/Users/pengyuyan/Desktop/singleFrame/images/singleImages/{}.png".format(i)color = cv2.imread(original_image_name)color = cv2.resize(color, (1214, 1354))width, height, ndim = color.shape[1], color.shape[0], color.shape[2]map_x = dest_array[:, 1].reshape((height, width)).astype(np.float32)map_y = dest_array[:, 0].reshape((height, width)).astype(np.float32)mapped_img = cv2.remap(color, map_y, map_x, cv2.INTER_LINEAR)cv2.imwrite("/Users/pengyuyan/Desktop/singleFrame/images/mytest_result_images/speed-optimization/INTER_LINEAR/INTER_LINEAR_result_{}.png".format(i), mapped_img)end = time.time()print("end!!! " + time.ctime())print("duration: {}".format(end - start))

最终的速度，4张图片的时间如下所示。

start!!! Fri Dec 13 11:41:56 2019
end!!! Fri Dec 13 11:41:56 2019
duration: 0.43956685066223145Process finished with exit code 0

一些有用的链接

https://www.youtube.com/watch?v=B7qrgrrHry0&feature=youtu.be
https://en.wikipedia.org/wiki/Distortion_(optics)

图像矫正：桶形畸变矫正的原理及python简易实现与加速相关推荐

python 取array并集_Python内置数据结构原理与性能简易分析
ins @ngladc 文末左下方阅读原文指向了本人博客链接,不含广告.参考资料中的相关链接,可以在博客文章的最下方获取.推荐苹果手机用户使用浅色模式观看. 前言对于一些算法题,可以使用Python ...
冲量（momentum）的原理与Python实现
冲量(momentum)的原理与Python实现前言参考:https://www.jianshu.com/p/58b3fe300ecb 梯度下降法(Gradient Descent)是机器学习中最 ...
梯度下降法快速教程 | 第三章：学习率衰减因子（decay）的原理与Python实现
北京 | 深度学习与人工智能 12月23-24日再设经典课程重温深度学习阅读全文> 正文共3017个字.11张图.预计阅读时间:8分钟前言梯度下降法(Gradient Descent)是 ...
梯度下降法快速教程 | 第二章：冲量（momentum）的原理与Python实现
北京 | 深度学习与人工智能研修 12月23-24日再设经典课程重温深度学习阅读全文> 01 前言梯度下降法(Gradient Descent)是机器学习中最常用的优化方法之一,常用来求解 ...
python gdbt+fm_GBDT回归的原理及Python实现
提到GBDT回归相信大家应该都不会觉得陌生(不陌生你点进来干嘛[捂脸]),本文就GBDT回归的基本原理进行讲解,并手把手.肩并肩地带您实现这一算法. 完整实现代码请参考本人的p...哦不是...git ...
手把手教你EMD算法原理与Python实现(更新)
Rose今天主要介绍一下EMD算法原理与Python实现.关于EMD算法之前介绍过<EMD算法之Hilbert-Huang Transform原理详解和案例分析>, SSVEP信号中含有自 ...
信号处理之频谱原理与python实现
目录频谱分析 FFT频谱分析原理下面就用python案例进行说明案例1 案例2 短时傅里叶变换STFT 本分享为脑机学习者Rose整理发表于公众号:脑机接口社区.QQ交流群:941473018 ...
倒频谱原理与python实现
目录倒频谱定义倒频谱python案例本教程为脑机学习者Rose发表于公众号:脑机接口社区 .QQ交流群:903290195 倒频谱定义倒频谱可以分析复杂频谱图上的周期结构,分离和提取在密集调频 ...
python实现逻辑回归的流程_逻辑回归原理及其python实现
September 28, 2018 7 min to read 逻辑回归原理及其python实现原理逻辑回归模型: $h_{\theta}(x)=\frac{1}{1+e^{-{\theta}^ ...

图像矫正：桶形畸变矫正的原理及python简易实现与加速