基于铅垂线的相机畸变标定算法

基于铅垂线的相机畸变标定算法

相机的径向畸变有两种方法来做标定，一种是标定板，一种是铅垂线法，本文主要是介绍关于铅垂线法来标定相机畸变。铅锤线法是指将相机画面中为弯曲的线段但是真实世界却是直线的线段通过一定的方法估计畸变方程的参数从而实现"拉"直。本文主要分为两小结，第一小结为介绍相机畸变产生的原因和畸变矫正模型，第二小结为基于铅垂线法相机畸变矫正的相关论文导读。

1. 相机畸变的原因和畸变矫正模型

相机畸变

相机的径向畸变是指光学系统对物体所成的像相对于物体本身而言是失真的，是光学透镜的固有特性，其原因是因为透镜边缘部分和中心部分的放大倍率不一致导致的。

径向畸变

许多广角镜头又很明显的径向畸变，对此不可能对他们做高精确度的图像重建。
摄像头的径向畸变(radial distortion) 分为三种类型: 桶形(barrel)、枕形（pincushion)以及鱼眼(fisheye)。

distortion_type.png

它导致一个圆形的畸变。
radial_distortion_result.png

偏心畸变

偏心畸变(decentering distortion) 是指相机镜头中心轴与成像平面不共线，导致图像中心位置的偏移，它通常出现在广角镜头中，因为广角镜头具有更大的视场角，需要更强的镜头弯曲来实现。

decentering.png

切向畸变

切向畸变(tangential distortion) 是指相机镜头和相机成像平面不完全平行，这回导致相机镜头的中心投影位置发生偏移，进而使得图像中的物体在图像边缘的位置发生形变产生切向畸变，这在网络相机上比较常见。

tangential_distortion.png

它会导致让图像径向畸变的圆形畸变转为椭圆型的畸变
tangential_distortion_result.png

畸变矫正模型

径向畸变矫正模型

解决相机径向畸变的方式有两种模型，一种是多项式乘法模型，一种是多项式除法模型。它们实现的原理就是以原画面中心点为圆心的不同半径上的圆上的像素点重新映射到新的一幅图片上不同的位置中。由于鱼眼的变化幅度很大，通常鱼眼相机的畸变矫正是单独采用另外的方式处理的。

多项式乘法模型
\begin{aligned} L(r) = 1 + k_1r^2+k_2r^4+...\end{aligned}
多项式除法模型
\begin{aligned}L(r)= \frac{1}{1 + k_1r^2+k_2r^4+...}\end{aligned}
其中:
\begin{aligned} x_c^{'} = x_c L(r) \\ y_c^{'} = y_cL(r)\\ r_c^2 = x_c^2 +y_c^2 \end{aligned}

上式中(x_c,y_c)表示的是原图片的像素位置，(x_c^{'},y_c^{'})表示的是畸变矫正吼的像素点，r_c表示的是(x_c,y_c)点到原图片中心点的距离。经过验证可知除法模型的有事在于用更小的阶数就能达到更高阶数多项式模型的矫正效果。

切向畸变矫正模型

如果相机有切向畸变的话则会添加一个切向畸变矫正模型。仅作了解，本文忽视这个畸变。

\begin{aligned} x_c^{'} = x+[2p_1xy+p_2(r^2+2x^2)] \\ y_c^{'} = y+[2p_2xy+p_1(r^2+2x^2)]\\ \end{aligned}

2. 畸变矫正文献摘要

本节内容前面3篇会介绍乘法模型在相机畸变矫正中的应用，后面2篇主要介绍的是除法模型在相机畸变的应用。这些畸变矫正算法的共性都是分为三个部分，第一部分是直线提取;第二部分是构造损失函数；第三部分是进行参数求解。有些文献的直线提取是手工提取，有的是用一些算法来取得弯曲的直线，参数求解部分以非线性解法较为流行。

2.1 Differential Methods for Nonmetric Calibration of Camera Lens Distortion

Moumen T. Ahmed and Aly A. Farag/University of Louisville/CVPR/2000

Motivation

摄像机的镜头畸变在中广角镜头上很明显，不同于之前的非度量(nonmetric)方法，本文针对畸变图片中直线的分析提出一种快速的封闭解(close-form solution)来得到畸变的参数。

Method

A. 曲线拉直

假设直线l是没有畸变图片上的直线，其上面的点为(x^u,y^u)，他们满足以下的方程
\begin{aligned}ax^u+by^u+c = 0\end{aligned}
其中a,b,c都是l的常量，并且l的系数为-\frac{a}{b}。在畸变图片中(x^d,y^d)也都在直线l上，因此上面的方程可以改写为:
\begin{aligned}f(x^d,y^d) = a x^u(x^d,y^d)+by^u(x^d,y^d)+c = 0\end{aligned}
对上式求导可得
\begin{aligned}\partial f = a[\frac{\partial x^u}{\partial x^d} \delta x^d+\frac{\partial x^u}{\partial y^d}\delta y^d]+b[\frac{\partial y^u}{\partial x^d} \delta x^d+\frac{\partial y^u}{\partial y^d}\delta y^d]=0\end{aligned}
又斜率-\frac{a}{b}要等于正切值\frac{\partial y^d}{\partial x^d}，联立可得在畸变图片中，直线l在(x^d,y^d)点的斜率s:
\begin{aligned}s(x^d,y^d) = \frac{\frac{\partial y^u}{\partial x^d} +\frac{\partial y^u}{\partial y^d}\frac{\delta y^d}{\delta x^d}}{\frac{\partial x^u}{\partial x^d} +\frac{\partial x^u}{\partial y^d}\frac{\delta y^d}{\delta x^d}}\end{aligned}
令(x^d_i,y^d_i),i=1,...,N是畸变图像曲线上的点，他们在矫正后得图像应属于同一直线，所以求解每个点间斜率得最小误差来求得畸变方程的参数
E = \sum_{i=2}^{N}(s(x^d_i,y^d_i)-s(x^d_{i-1},y^d_{i-1}))^2

B. 畸变标定

上面的损失方程可以使用非线性(non-linear)的解法去求最终的结果来得到更准确的结果，但是它会花费更长的时间以及如果遇到一个比较差的初始值可能会得到错误的结果或是局部最优。封闭解(closed-form)会比非线性解更精确且速度更快。

a. 求根解法
当畸变矫正方程的未知变量p大于2的时候没有解析解;当p中只有一个K1和P1的未知数的时候有封闭解。
在求解过程中一般将点分割成M个groups，计算每个group的畸变系数，然后针对这一系列的根R_m = {r_{mj} | j = 1,...,n_m}然后寻找到一个根r使得所有根能最小化
\begin{aligned}\sum_{m=1}^M d(r,R_m) = \sum_{m=1}^M min |r-r_{mj}| \end{aligned}

b. 线性解法
当p<=2的时候封闭解，形如AX=B,A是一个N.L * p的矩阵，L是曲线数量，N则是每段线段用来计算的点数量。X是未知的畸变参数p未知参数的数量。通过奇异矩阵分解就可以得到特征根。

c. 畸变中心的不确定性
线性解需要确定好畸变中心,所以作者提出两个引理，引理的具体证明在附录上
引理1: 畸变中心在图像中心处(c_x,c_y)偏移了(\Delta_x,\Delta_y)等于在式子中添加了中心偏移的式子-K_1\Delta_x和-K1\Delta_y
引理2:在镜头偏心畸变的影响下，(c_x,c_y)偏移了(\Delta_x,\Delta_y)对系数的影响不大。

Experiments

实验的图片大小为320x242。


DMFN_Syn.png

由上可知，非线性的解法得到的参数会更加接近合成图片的畸变参数。


DMFN_real.png

在真实图片中，可见是都把弯曲的线段拉直了。

Conclusion

本篇文章提出两种标定相机畸变参数的方法。在古早时代，非线性解在SGI-02的机器上需要运行3.2s但是封闭解只需要0.05s，但是这点算力在今天的机器上不差那么点。

2.2 Nonmetric Lens Distortion Calibration: Closed-form Solutions, Robust Estimation and Model Selection

Moumen T. Ahmed and Aly A. Farag/University of Louisville/ICCV/2003

Motivate

之前的相机标定方法需要使用者参与标定的过程，本文提出一种基于最小中值二乘法(the-least-median-of-squares (LMedS) estimator)的自动标定方法，基于此方法有三个优点:
(1). 通过LMedS的方法，不需要人工介入选取真实世界是直线但是畸变图片上曲线的线段。
(2). 通过线性的方法求得closed-form solution后作为非线性解的初始值，能更好解决畸变参数卡在局部最优值的问题。
(3). 通过几何推断来选取较优的模型。

Method

参数优化

本文的参数优化是[1]中的非线性解和线性解的结合，具体可参考[1]。

边缘点的筛选

要实现自动标定需要提取图片上的曲线，但是提取图片上的曲线会遇到两个问题，一、提取的曲线在现实世界不是直线；二、一条长曲线断成多截的线段导致有用信息偏少。基于上述的原因本文通过最小中值平方法来去除错误的外点。

算法流程

(1). 提取边缘线段，连成直线并用最小中值平方法过滤点。
(2). 使用bucketing techniques来将图片分成b*b的区域，在每个bucket区域采样一半的点。
(3). 多次采样计算每次的畸变模型参数。
(4). 通过Akaike’s AIC 和 Rissanen’s MDL方法来选取其中比较好的畸变矫正模型。

Experiments

不同噪声情况下下算法矫正后几何推断的指标。


NLDC_geometric.png

论文示例图线段提取算法以及矫正的效果都看着还不错，c图是封闭解的优化结果，d图是与非线性解的结果。


NLDC_real.png

Conclusion

对曲线拉直的畸变矫正方法没什么改进，只是基于作者在2000的ICCV文章上增加了一些提取直线和增加模型矫正效果的评判指标，意义不是很大。

2.3 Line-Based Correction of Radial Lens Distortion

B. PRESCOTT AND G. F. MCLEAN/University of Victoria/GRAPHICAL MODELS AND IMAGE PROCESSING/1995

Motivation

在此之前对相机的标定需要使用标定板来标定许多参数(内参、外参、畸变参数)来解决径向畸变的问题。作者提出了一种不需要标定板就可以做好畸变矫正的方法。
相机的标定可以分为三个部分，1. 外参(相机的朝向、位置等);2. 内参(焦距等);3. 畸变参数(径向畸变和偏心畸变)，本文在不适用标定板的情况下只对相机的畸变参数进行标定。

Method

直线提取算法

作者使用了[4]中的Line Support Region(LSR)直线提取算法，该算法对局部梯度方向的边缘进行了连接，从而得到真实世界的直线在畸变图片中为曲线的线段。

参数优化与更新

令(x,y)为真实世界中的直线上的点，满足
\begin{aligned}x \cos \theta +y \sin \theta = \rho\end{aligned}
令LSR_m = \{(x^{'}_{k},y^{'}_{k}),k=1,...,N_m\}直线提取算法提取的第m条线段,相应的点经过畸变矫正后为(x^u_k,y^u_k)满足
\begin{aligned}x^u_k\cos\theta_m+y^u_k\sin\theta_m \approx \rho_m\end{aligned}
对于LSR算法提取的直线都能满足上面的方程时，此时的畸变参数就是最优解
\begin{aligned}F = \sum^M_{m=1}\sum^{N_m}_{k-1}(x^u_k \cos \theta_m +y^u_k \sin \theta_m -\rho_m)^2\end{aligned}
作者将K_1,K_2的初始值都设置为0，c_x,c_y设置为图像的中心点。

Experiments

由于畸变矫正的评价指标不明，作者使用了合成和真是的图片来测。为了体现算法的有效性，作者定义了一个曲面指标，计算原始图片上像素点与欧氏距离。
\xi_d (x,y)= \sqrt{(x-x^{'})^2+(y-y{'})^2}
以下是作者合成的图片与经过算法畸变矫正后的图片

LBCO_Syn.png

未矫正和矫正后的指标可视化的结果如下:
LBCO_surface.png

由上可知，经过矫正后的图片在改指标上有很大的改进。

Conclusion

本文的实验图片偏简单一些，实验的图片也比较少。由于代码没开源，具体的计算细节无从得知。

2.4 Automatic Lens Distortion Correction Using One-Parameter Division Models

Miguel Alem´an-Flores1/Universidad de Las Palmas de Gran Canaria/IPOL/2013 17页好多

Introduction

作者使用了添加了畸变参数的霍夫变换方法来更好抽取直线，然后解一项式的除法模型参数来做畸变矫正。
a. 使用Canny来做边缘检测
b. 通过使用改进的霍夫变换方法来得到畸变的初始值以及投票数最多的直线。
c. 畸变参数优化
d. 图像畸变矫正

Method

A. 使用canny算法获取边缘点以及边缘点的角度

就是原本的canny算法改了一些超参。

B. 改进式霍夫变换(Improved Hough Transform)

a. 畸变参数变换为分辨率无关
畸变方程L(r)中的k_1的值和r相关，以至于当图片的大小变化的时候k_1的值也要跟着变换。
作者提出一个新的参数p=(r^{'}_{max} - r_{max})/r_{max}，其中r^{'}_{max}为畸变校正的距离中心点最远端的值，r_{max}是畸变图像中最远端的值。因此，将上式子带入到畸变方程中，可得:
\begin{aligned} k_1 = \frac{-p}{(1+p)r^2_{max}}\end{aligned}
对畸变方程求r导倒数，由于随着r的增加r^{'}也会增加，可知畸变方程是一个单调递增的函数。
\frac{(1+k_1r^2)-2k_1r^2}{(1+k_1r^2)^2} >0 \Rightarrow k_1<\frac{1}{r^2_{max}}
因此:
p>-0.5

b. 3D霍夫变换
(1). 对于上式中得p，作者选取一个间隔[p_{min},p_{max}]其中(p_{min})>-0.5 对于区间内得所有p值可表示为:
p \in P = \{ p_{min}+n\delta_p : n = 0,1,...,\frac{p_{max}-p_{min}}{\delta_p}\}，\delta_p = 0.1
因此作者将霍夫变换的空间变换成(d,\alpha,p)三维空间。
(2). 针对不同的p值作者会重新计算边缘点的位置以及方向。
(3). 对于变换后边缘点要与其相符的直线角度小于\delta_{\alpha}度。
(4). 对于所有满足c步条件的点，求得平均距离d_l,剔除与d_l与绝对值相差2的点，设置投票的权重设置为\frac{1}{1+d_p}
最终p与投票分数的关系如图

ALDC_p.png

有图可知，p的值比较小的时候能有比较多的点在直线上。

C. 畸变参数优化

a. 作者移除小于5个点的直线，并且选取最长的Nl条线段。因为小直线不能提供比较有用的相关信息，并且更大的直线更有可能在真实世界中是直线，选取的直线数量是30条。
b. 对于p \in R 和集合中的直线的点\{(x_{ji},y_{ji})\}和对应的经过矫正的点\{x^{'p}_{ji},y^{'p}_{ji}\}，每条直线的损失函数如下
D(\alpha,d) = \sum_{i=1}^{N(j)}(\cos(\alpha)x^{'p}_{ji} + \sin(\alpha)y^{'p}_{ji}+d)^2
对于所有的直线得到的损失方程为:
E(p) = \frac{\sum^{Nl}_{j=1}\sum^{N(j)}_{i=1}(\cos(\alpha^p_j)x^{'p}_{ji}+\sin(\alpha^p_j)y^{'p}_{ji}+d^p_j)^2}{\sum^{Nl}_{j=1}N(j)}
其中p的初始值为p0，求导更新优化。

Experiments

作者演示了两张图片，一张是标定板，一张是建筑。由下图可知，作者的直线提取算法还是挺不错的效果，图上的指向基本上都被拉直了。


ALDC_real.png

Conclusion

比较遗憾的是也是没有一个指标来表明该算法的效果以及和同类型的算法进行对比。

2.5 An Iterative Optimization Algorithm for Lens Distortion Correction Using Two-Parameter Models

Miguel Alem´an-Flores1/Universidad de Las Palmas de Gran Canaria/IPOL/2016 39页好多

Introduction

本文是[5]的作者出的改进版。作者在本文的工作中做了以下的改进:
a. 增加了求解的畸变模型参数。优化了[5]中的基于除法模型的一个参数的自动畸变矫正算法为两个畸变参数和畸变中心点。
b. 对canny算法提出来的边缘点做了筛选，去除了一部分非直线的边缘点。
c. 作者对改进式霍夫变换提取出来的直线做了后处理，让一些小线段可以连接起来。
d. 作者通过增加一些边缘点的策略多次迭代求解最优参数。
e. 图像畸变矫正

Method

A. 使用Canny算法做边缘检测

使用canny算法改些超参数得到边缘点。对提取得边缘点做进一步处理来尽量去除非直线的边缘点。
作者通过定义一个EOSM(p)的指标来衡量边缘点与其周边相邻的边缘点方向上差异，两点间角度差异越大\cos(\alpha_p - \alpha_q)的值越小
\begin{aligned}EOSM(p) = \sum_{q \in Edge points \bigcap Neighborhood(p)}\cos(\alpha_p)cos(\alpha_q)+\sin(\alpha_p)\sin(\alpha_q)\\ = \sum_{q \in Edge points \bigcap Neighborhood(p)} cos(\alpha_p-\alpha_q)\end{aligned}
首先移除EOSM(p)小于一定阈值的点，之后为了加速后面的步骤的运算量，作者去除邻域中EOSM(p)最小的点来减少点的数量。
(插入图片)

B. 3D霍夫变换

前一部分和[5]的步骤相同。作者在这片文章中增加一个后处理的步骤。闭合相邻的小线段以及最后去除孤立很小的线段。这部分繁琐又复杂，具体细节得看看代码。

C. 迭代式参数更新

作者将第一个参数的初始值设置为霍夫变换得到得最佳参数。
a. 参数更新
令d为要优化的参数(k_1,k_2,x_c,y_c)，Nl为直线的数量，N(j)表示为第j条直线的点，C_d(x^{'}_{ji})表示为经过畸变矫正模型后的的第j条线，则要优化的目标(能量)函数为:
E(d) = \sum^{Nl}_j\sum^{N(j)}_i distance(C_d(x^{'}_{ji}),line_j)^2
E(d)在d_0处的泰勒展开式为:
E(d) = E(d_0) +(\nabla E(d_0))^T(d-d_0)+\frac{1}{2}(d-d_0)^T (\nabla)^2E(d_0)(d-d_0)+....
因为要最小值这个损失函数，所以前面泰勒展开式部分要等于0
\nabla E(d) \approx \nabla E(d_0)+(\nabla)^2E(d_0)(d-d_0)=0
通过Gauss method可求得
d_{n+1} = d_n+((\nabla)E(d_n)+\gamma Id)^{-1}(-\nabla E(d_n))
其中\gamma是为了保证随着参数的优化减小(学习率的作用).最终k_1和k_2与p1,p2的关系如下
k_1 = \frac{\frac{-p_1}{1+p_1}+\frac{16p_2}{1+p_2}}{-12r^2_2},k_2 = \frac{\frac{-4p_2}{1+p_2}+\frac{p_1}{1+p_1}}{-12r^4_2}

完整的参数求解过程详见作者的另外一篇论文[8]
b. 迭代
在参数优化的过程中，作者还将霍夫变换后处理中剔除的点经过新参数校验后加入计算的点集中，这个做法可以检测出更多的直线以及直线上又更多的点在畸变很严重的时候又很好的效果。

Expirements

由于没有一个比较指标，实验表明在图片上，两个参数的除法模型(Div2pIO)比一个参数的除法模型(Div1p)掉落在直线上的点会更多。作者的实验图片，相对于之前的那一篇使用了更大的焦距使得图片形变更大，在d图上，可以知道，提取的直线也基本上被拉直了。


AIOA_scene.png

AIOA_scene_table.PNG

Conclusion

这篇文章在判断点是否在直线上下了很大的功夫。作者的实验也表明了，同样参数的除法模型畸变纠正的效果比乘法模型好。

参考文献

[1] Differential Methods for Nonmetric Calibration of Camera Lens Distortion
[2] Nonmetric Lens Distortion Calibration: Closed-form Solutions, Robust Estimation and Model Selection
[3] Line-Based Correction of Radial Lens Distortion
[4] Extracting straight lines
[5] Automatic Lens Distortion Correction Using One-Parameter Division Models
[6] An Iterative Optimization Algorithm for Lens Distortion Correction Using Two-Parameter Models
[7] Computer Vision: Algorithms and Applications 2nd Edition
[8] Invertibility and Estimation of Two-Parameter Polynomial and Division Lens Distortion Models

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 201,468评论 5赞 473
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 84,620评论 2赞 377
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 148,427评论 0赞 334
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,160评论 1赞 272
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,197评论 5赞 363
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,334评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 37,775评论 3赞 393
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,444评论 0赞 256
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 40,628评论 1赞 295
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,459评论 2赞 317
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,508评论 1赞 329
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,210评论 3赞 318
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 38,767评论 3赞 303
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,850评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,076评论 1赞 258
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 42,627评论 2赞 348
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,196评论 2赞 341

推荐阅读更多精彩内容