CN103258327A - 一种基于二自由度摄像机的单点标定方法 - Google Patents

Abstract

本发明属于机器视觉领域，公开了一种基于二自由度摄像机的单点标定方法，包括根据标识点在基准坐标系中的三维物理坐标和与标识点对应的像点在图像上的二维像素坐标获得标定的关系模型；采用线性方法计算标定的关系模型获得水平转角和垂直转角。本发明提供的这种标定方法在求解过程中，巧妙的将复杂费时的非线性方程组的求解问题转化为带参数的线性方程组的求解问题，并且在各种情形下给出了水平转角的正、余弦闭式解和垂直转角的正、余弦闭式解；这种单点标定方法不仅是在只有一个标识点，且需要在线标定的应用环境中的唯一选择。且通过实验表明，采用本发明提供的这种单点标定方法得到的PTZ摄像机的水平转角Pan和垂直转角Tilt的绝对误差也非常低。

Description

一种基于二自由度摄像机的单点标定方法

技术领域

本发明属于机器视觉领域，更具体地，涉及一种基于二自由度摄像机的单点标定方法。

背景技术

目前，Pan-Tilt-Zoom(PTZ)摄像机在越来越多的场合下被应用于视频采集系统。PTZ摄像机可以实现镜头的水平和垂直方向的转动，可以调整摄像机的焦距。在采用定焦摄像机的应用环境中，PTZ摄像机只有水平方向上的转角Pan和垂直方向上的转角Tilt两个自由度。若利用这种视频采集系统在线进行三维立体测量，就需要对摄像机进行在线标定，以标定出PTZ摄像机的水平方向上的转角Pan和垂直方向上的转角Tilt。标定结果的准确性直接影响了三维立体测量的精度。

在摄像机的使用过程中，若摄像机的内参数保持不变，则可以在实验室条件下预先标定摄像机的内参数。

若摄像机的位置设定之后，不再进行调整，则摄像机在基准坐标系中的三维物理坐标也可以通过各种方式事先确定。

低精度PTZ摄像机可以输出精度为0.5°的摄像机的转角Pan和Tilt的整数值角度，其中，Pan的角度范围在0°～360°之间，Tilt的角度范围在-78°～0°之间。显然，这种整数值精度的Pan和Tilt值不能满足三维精确测量的需求。

为了给出PTZ摄像机的转角Pan和Tilt的精确值，就需要采用传统标定方法、自标定方法和基于主动视觉的标定方法等常用的摄像机标定方法来标定Pan和Tilt。

传统的摄像机标定方法是利用标定块与拍摄图像之间的约束关系建立关于摄像机待标定参数的标定模型，通过优化算法来求解摄像机的待标定参数。传统的标定方法标定结果准确度高，但是所采用的优化算法费时，需要较多的标识点，并且标识点所在的标定块一般比较昂贵。

自标定方法不需要特定的标定块，而是利用环境中的刚体型，通过对比多幅图像中的刚体型上的特定几何结构来标定摄像机的待标定参数。自标定方法不需要标定块，但是在缺少刚体型的应用环境中会失效，而且算法的鲁棒性和稳定性都差。

主动视觉的标定方法是将摄像机精确安装于精密可控平台，在已知摄像机的运动信息基础上，利用图像上同名点的对应关系对摄像机的待标定参数进行标定。该方法鲁棒性强，但无法应用在摄像机运动未知的场合，所需平台精度要求较高，成本也较高。

当应用环境中没有特定几何结构的刚体型，无法获取摄像机的精确运动信息，也没有标定块时，上述三种标定方法均不适用。

张正友提出的标定方法必须利用标定板上至少六个标识点才能进行摄像机的标定，可在实验室环境中用来对摄像机内参数进行标定。但不适用于在室外环境中对摄像机外参数的在线标定。

Jacek Komorowski和Przemyslaw Rokita提出的外参数标定方法，需要在不同的未知姿态下获取一系列图像。根据多幅图像中多个同名点的对应关系，采用迭代的方法获得这些不同姿态之间的相对姿态信息。由于采用了复杂并且费时的迭代方法，这种标定方法不适用于在线标定。

在户外恶劣的陌生环境中，有时无法获得足够数目的标识点的三维物理坐标信息。极端情形下，只有一个标识点的三维物理坐标信息可被用于Pan和Tilt的标定。在这种仅有一个标识点且需要在线定标的应用环境中，上述现有的像机标定方法都不适用。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种基于二自由度摄像机的单点标定方法，其目的在于仅利用唯一的单个标识点可准确的标定出PTZ摄像机的水平转角和垂直转角，由此解决现有技术中针对只有一个标识点无法对摄像机标定的技术问题。

为实现上述目的，按照本发明的一个方面，提供了一种基于二自由度摄像机的单点标定方法，包括下述步骤：

S1：根据标识点在基准坐标系中的三维物理坐标和与所述标识点对应的像点在图像上的二维像素坐标获得标定的关系模型；

S2：采用线性方法计算所述标定的关系模型获得水平转角和垂直转角。

具体地，所述标定的关系模型为：

$Z\left(\begin{array}{c}U\\ 1\\ -V\end{array}\right)=\left(\begin{array}{ccc}\cos \left(\mathrm{\ΔP}\right)& \sin \left(\mathrm{\ΔP}\right)& 0\\ -\sin \left(\mathrm{\ΔP}\right)\cos \left(\mathrm{Tilt}\right)& \cos \left(\mathrm{\ΔP}\right)\cos \left(\mathrm{Tilt}\right)& \sin \left(\mathrm{Tilt}\right)\\ \sin \left(\mathrm{\ΔP}\right)\sin \left(\mathrm{Tilt}\right)& -\cos \left(\mathrm{\ΔP}\right)\sin \left(\mathrm{Tilt}\right)& \cos \left(\mathrm{Tilt}\right)\end{array}\right)\left(\begin{array}{c}{X}_p\\ Y_p\\ Z_p\end{array}\right);$

U表示标识点在焦距为1时的正则化像平面上所对应像点的像素列坐标，V表示标识点在焦距为1时的正则化像平面上所对应像点的像素行坐标，Z表示标识点P在像空坐标系C-XYZ中的第三维坐标，(X_p，Y_p，Z_p)^T表示标识点在临时坐标系C-x’y’z中的三维物理坐标；ΔP表示水平转角与初始参考值之间的差值，Tilt表示垂直转角。

具体地，所述临时坐标系C-x’y’z的获取方法包括：将基准坐标系O-xyz平移得到过渡临时坐标系C-xyz；再将过渡临时坐标系C-xyz绕其Cz轴正向右手螺旋旋转角度P₀后得到所述临时坐标系C-x’y’z。

具体地，所述步骤S2包括：

S21：将所述标定的关系模型转化成带水平转角参数并关于垂直转角的正、余弦的线性关系模型；

S22：根据所述线性的关系模型获得垂直转角的正、余弦闭式表达式；

S23：根据垂直转角的正、余弦闭式表达式并利用三角函数约束得到水平转角的正切的闭式表达式；

S24：根据所述水平转角的正切闭式表达式和所述垂直转角的正、余弦闭式表达式得到水平转角和垂直转角。

具体地，所述线性关系模型为：

$\left(\begin{array}{cc}U{Y}_p\cos \left(\mathrm{\ΔP}\right)-U{X}_p\sin \left(\mathrm{\ΔP}\right)& Z_{p}U\\ V{X}_p\sin \left(\mathrm{\ΔP}\right)-V{Y}_p\cos \left(\mathrm{\ΔP}\right)-Z_p& Y_p\cos \left(\mathrm{\ΔP}\right)-V{Z}_p-X_p\sin \left(\mathrm{\ΔP}\right)\end{array}\right)\left(\begin{array}{c}\cos \left(\mathrm{Tilt}\right)\\ \sin \left(\mathrm{Tilt}\right)\end{array}\right)$

$=\left(\begin{array}{c}{X}_p\cos \left(\mathrm{\ΔP}\right)+Y_p\sin \left(\mathrm{\ΔP}\right)\\ 0\end{array}\right).$

具体地，当线性关系模型的系数矩阵行列式不等于零时，垂直转角的正、余弦闭式表达式为：

$\cos \left(\mathrm{Tilt}\right)=\frac{(Y_p\cos \left(\mathrm{\ΔP}\right)-V{Z}_p-X_p\sin \left(\mathrm{\ΔP}\right))(X_p\cos \left(\mathrm{\ΔP}\right)+Y_p\sin \left(\mathrm{\ΔP}\right))}{U[{Z_p}^2+{(Y_p\cos \left(\mathrm{\ΔP}\right)-X_p\sin \left(\mathrm{\ΔP}\right))}^2]}$ 和

$\sin \left(\mathrm{Tilt}\right)=\frac{(Z_p+V{Y}_p\cos \left(\mathrm{\ΔP}\right)-V{X}_p\sin \left(\mathrm{\ΔP}\right))(X_p\cos \left(\mathrm{\ΔP}\right)+Y_p\sin \left(\mathrm{\ΔP}\right))}{U[{Z_p}^2+{(Y_p\cos \left(\mathrm{\ΔP}\right)-X_p\sin \left(\mathrm{\Δp}\right))}^2]};$

水平转角的正切闭式表达式为：

$\tan \left(\mathrm{\ΔP}\right)=\frac{-X_p{Y}_p(U^2+V^2+1)\text{\±U}\sqrt{(X_p^2+Y_p^2+Z_p^2)[(X_p^2+Y_p^2)(V^2+1)-U^2{Z}_p^2]}}{(V^2+1)Y_p^2-U^2(X_p^2+Z_p^2)}$ 或

$\tan \left(\mathrm{\ΔP}\right)=\frac{U^2{Z}_p^2+U^2{Y}_p^2-(V^2+1)X_p^2}{2X_p{Y}_p(U^2+V^2+1)}.$

具体地，当线性关系模型的系数矩阵行列式等于零时，所述垂直转角的正切闭式表达式为 $\tan \left(\mathrm{Tilt}\right)=\frac{Z_p+V{Y}_p\cos \left(\mathrm{\ΔP}\right)-V{X}_p\sin \left(\mathrm{\ΔP}\right)}{Y_p\cos \left(\mathrm{\ΔP}\right)-X_p\sin \left(\mathrm{\ΔP}\right)-V{Z}_p};$ 所述水平转角的正切闭式表达式为 $\tan \left(\mathrm{\ΔP}\right)=-\frac{X_p}{Y_p}.$

总体而言，本发明提供的这种标定方法在求解过程中，巧妙的将复杂费时的非线性方程组的求解问题转化为带参数的线性方程组的求解问题，并且在各种情形下给出了sin(ΔP)、cos(ΔP)、sin(Tilt)和cos(Tilt)的闭式解。这种单点标定方法不仅是在只有一个标识点，且需要在线标定的应用环境中的唯一选择。且通过实验表明，采用本发明提供的这种单点标定方法得到的PTZ摄像机的水平转角Pan和垂直转角Tilt的绝对误差也非常低。

附图说明

图1是本发明实施例提供的基准坐标系O-xyz、二维像素坐标系和像空坐标系C-XYZ的示意图；

图2是本发明实施例提供的坐标系旋转示意图；

图3是本发明实施例提供的基于二自由度摄像机的单点标定方法实现流程图；

图4是本发明实施例提供的基于二自由度摄像机的单点标定方法中步骤S2的子流程图；

图5是本发明实施例提供的水平转角Pan角的绝对误差结果图；

图6是本发明实施例提供的垂直转角Tilt角的绝对误差结果图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明提供的基于二自由度摄像机的单点标定方法可以广泛应用于机器人导航系统、三维重建技术、生物医学领域以及3D虚拟技术中。

在摄像机的内参数矩阵、摄像机在基准坐标系中的三维物理坐标、一个标识点在基准坐标系中的三维物理坐标已知，摄像机只有水平和垂直方向的转动，并且需要在线标定出摄像机的水平转角Pan和垂直转角Tilt的应用环境下，我们提出了一种标定摄像机的外参数的方法，称之为“PTZ摄像机的单点在线标定方法”。摄像机标定是一种求解物点的三维物理坐标与其所对应像点的二维图像坐标之间的关系模型中的待定模型参数的过程。

如图1所示，PTZ摄像机的像空坐标系C-XYZ的原点C为摄像机的光心，坐标轴CX的正向与图像的列增加的方向一致，坐标轴CY的正向与图像的行增加的方向一致，坐标轴CZ的正向由摄像机光心指向图像平面Auv，并且与摄像机的光轴方向一致。

如图2所示，基准坐标系O-xyz和像空坐标系C-XYZ之间的变换关系可以描述为：(1)将基准坐标系0-xyz平移得到临时坐标系C-xyz，使得临时坐标系的坐标系原点为点C；(2)将临时坐标系C-xyz绕Cz轴正向右手螺旋转动Pan角，得到临时坐标系C-Xy′z；(3)将临时坐标系C-Xy′z绕CX轴正向右手螺旋转动(Tilt-90度)角，得到PTZ摄像机的像空坐标系C-XYZ。根据小孔成像原理，物点P在基准坐标系O-xyz中的三维物理坐标(x，y，z)^T，与其所对应的像点Q在图像上的二维像素坐标(u，v)^T之间的关系模型为：

$Z\left(\begin{array}{c}u\\ v\\ 1\end{array}\right)=\mathrm{KR}\left(\begin{array}{c}x-x_c\\ y-y_c\\ z-z_c\end{array}\right);---\left(1\right)$

其中，K是摄像机的内参数矩阵，R是基准坐标系0-xyz与PTZ摄像机的像空坐标系C-XYZ之间的旋转变换矩阵，(x_c，y_c，z_c)^T是摄像机在基准坐标系的三维物理坐标，Z是物点P在坐标系C-XYZ中的第三维坐标。

由图2所示的基准坐标系0-xyz和像空坐标系C-XYZ之间的变换关系知，R可表示为：

$R=\left(\begin{array}{ccc}1& 0& 0\\ 0& \sin \left(\mathrm{Tilt}\right)& -\cos \left(\mathrm{Tilt}\right)\\ 0& \cos \left(\mathrm{Tilt}\right)& \sin \left(\mathrm{Tilt}\right)\end{array}\right)\left(\begin{array}{ccc}\cos \left(\mathrm{Pan}\right)& \sin \left(\mathrm{Pan}\right)& 0\\ -\sin \left(\mathrm{Pan}\right)& \cos \left(\mathrm{Pan}\right)& 0\\ 0& 0& 1\end{array}\right)---\left(2\right)$

将(2)式代入(1)式，可得：

$Z\left(\begin{array}{c}u\\ v\\ 1\end{array}\right)=K\left(\begin{array}{ccc}1& 0& 0\\ 0& \sin \left(\mathrm{Tilt}\right)& -\cos \left(\mathrm{Tilt}\right)\\ 0& \cos \left(\mathrm{Tilt}\right)& \sin \left(\mathrm{Tilt}\right)\end{array}\right)\left(\begin{array}{ccc}\cos \left(\mathrm{Pan}\right)& \sin \left(\mathrm{Pan}\right)& 0\\ -\sin \left(\mathrm{Pan}\right)& \cos \left(\mathrm{Pan}\right)& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}x-x_c\\ y-y_c\\ z-z_c\end{array}\right)---\left(3\right)$

其中，(3)式给出了物点P在基准坐标系O-xyz中的三维物理坐标(x，y，z)^T，与其所对应的像点Q在图像上的二维像素坐标(u，v)^T之间的关系模型。其中，PTZ摄像机的内参数矩阵K，PTZ摄像机在基准坐标系中的三维物理坐标(x_c，y_c，z_c)^T，以及PTZ摄像机的水平转角Pan和垂直转角Tilt是模型参数。

在PTZ摄像机的内参数矩阵K，和PTZ摄像机在基准坐标系中的三维物理坐标(x_c，y_c，z_c)^T都已知的情形下，利用唯一的标识点在基准坐标系中的已知的三维物理坐标(x，y，z)^T，和该标识点在图像上所对应像点的已知的二维像素坐标(u，v)^T，来在线求解(3)式中的模型参数Pan和Tilt。

(1)重新构建模型：

根据标识点P在基准坐标系O-xyz中的三维物理坐标(x，y，z)^T和与其所对应的像点Q在图像上的二维像素坐标(u，v)^T获得标定的关系模型；

当PTZ摄像机的垂直转角Tilt变化范围在[-78°，0°]上时，tan(Tilt)一定为负，只要确定了tan(Tilt)，就能根据如下公式唯一的确定sin(Tilt)和cos(Tilt)；即： $\cos \left(\mathrm{Tilt}\right)=\frac{1}{\sqrt{1+{\tan }^2\left(\mathrm{Tilt}\right)}}$ 和 $\sin \left(\mathrm{Tilt}\right)=\frac{\tan \left(\mathrm{Tilt}\right)}{\sqrt{1+{\tan }^2\left(\mathrm{Tilt}\right)}}.$ 而水平转角Pan变化范围在[0°，360°]上，即使确定了tan(Pan)，也不能唯一的确定sin(Pan)和cos(Pan)。因此，为减少不确定性，我们以PTZ摄像机输出的整数形式的Pan角度P₀作为精确解算Pan的初始参考值，将P₀和精确的水平转角Pan的差值记为ΔP，有：Pan＝P₀+ΔP。其中，ΔP变化范围在[-0.5°，0.5°]上。这样，只要确定了tan(ΔP)，就能根据如下公式唯一的确定sin(ΔP)和cos(ΔP)；即： $\cos \left(\mathrm{\ΔP}\right)=\frac{1}{\sqrt{1+{\tan }^2\left(\mathrm{\ΔP}\right)}},$ $\sin \left(\mathrm{\ΔP}\right)=\frac{\tan \left(\mathrm{\ΔP}\right)}{\sqrt{1+{\tan }^2\left(\mathrm{\ΔP}\right)}}.$ 为了简化表达，记

$\left(\begin{array}{c}U\\ 1\\ -V\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 1\\ 0& -1& 0\end{array}\right)K^{-1}\left(\begin{array}{c}u\\ v\\ 1\end{array}\right),$ $\left(\begin{array}{c}{X}_p\\ Y_p\\ Z_p\end{array}\right)=\left(\begin{array}{ccc}\cos \left(P_0\right)& \sin \left(P_0\right)& 0\\ -\sin \left(P_0\right)& \cos \left(P_0\right)& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}x-x_c\\ y-y_c\\ z-z_c\end{array}\right)---\left(4\right);$

其中，(U，V)^T是标识点在焦距为1时的一种正则化像平面上所对应像点的像素坐标，(X_p，Y_p，Z_p)^T表示标识点在一个临时坐标系C-x’y’z中的三维物理坐标。将基准坐标系O-xyz平移为临时坐标系C-xyz后，再将临时坐标系C-xyz绕Cz轴正向右手螺旋旋转角度P₀后，就会形成这个临时坐标系C-x’y’z。

这样，式(3)可以改写为：

$Z\left(\begin{array}{c}U\\ 1\\ -V\end{array}\right)=\left(\begin{array}{ccc}\cos \left(\mathrm{\ΔP}\right)& \sin \left(\mathrm{\ΔP}\right)& 0\\ -\sin \left(\mathrm{\ΔP}\right)\cos \left(\mathrm{Tilt}\right)& \cos \left(\mathrm{\ΔP}\right)\cos \left(\mathrm{Tilt}\right)& \sin \left(\mathrm{Tilt}\right)\\ \sin \left(\mathrm{\ΔP}\right)\sin \left(\mathrm{Tilt}\right)& -\cos \left(\mathrm{\ΔP}\right)\sin \left(\mathrm{Tilt}\right)& \cos \left(\mathrm{Tilt}\right)\end{array}\right)\left(\begin{array}{c}{X}_p\\ Y_p\\ Z_p\end{array}\right)---\left(5\right)$

PTZ摄像机的单点在线定标的目标就是，能够利用式(5)中已知的(X_p，Y_p，Z_p)^T和(U，V)^T，消去未知参数Z，解算出待标定参数ΔP和Tilt。其中，ΔP属于[-0.5°，0.5°]，Tilt属于[-78°，0°]。

(2)将非线性方程组求解问题转化为线性方程组的求解问题，即采用线性方法计算上述标定的关系模型得到水平转角和垂直转角。

(2.1)将上述标定的关系模型转化成带水平转角参数ΔP并关于垂直转角Tilt的正、余弦的线性关系模型；

消去式(5)中的未知参数Z，我们可以得到：

U(-X_psin(ΔP)cos(Tilt)+Y_pcos(ΔP)cos(Tilt)+Z_psin(Tilt))

＝X_pcos(ΔP)+Y_psin(ΔP)

-V(-X_psin(ΔP)cos(Tilt)+Y_pcos(ΔP)cos(Tilt)+Z_psin(Tilt))

＝X_psin(ΔP)sin(Tilt)Y_pcos(ΔP)sin(Tilt)+Z_pcos(Tilt)                (6)

式(6)中，U，V，X_p，Y_p，Z_p都是已知的标量，ΔP和Tilt是待标定的未知参数。显然，式(6)是由2个关于未知参数ΔP和Tilt的三角函数型非线性方程构成的非线性方程组。

已知一个标识点的信息就能够列出两个如(6)式所示的方程。由形如(6)式的非线性方程组最多可以求解出其中的两个未知参数。理论上讲，利用唯一的标识点，构建出式(6)后，就可以采用非线性方程组的数值解法求解出未知参数ΔP和Tilt。但是，非线性方程组可能在求解的过程中出现无解或者多解的情况，并且用于求解非线性方程组的数值解法往往都需要多次迭代，比较复杂费时，不能满足在线定标的要求。因此，必须寻找更高效的求解方法。

观察式(6)，可以注意到，若将sin(Tilt)、cos(Tilt)、sin(ΔP)、cos(ΔP)分别看作4个变量，则式(6)可以看作是由2个这4个变量的2次方程构成的关于这4个变量的2次方程组。显然，sin(Tilt)、cos(Tilt)、sin(ΔP)、cos(ΔP)还必须要满足平方和为1的2个约束。这样，可将由式(6)求解ΔP和Tilt的非线性方程组的求解问题，转化为求解4个关于4个变量sin(Tilt)、cos(Tilt)、sin(ΔP)、cos(ΔP)的2次方程的2次方程组的求解问题。相比于直接求解关于未知参数ΔP和Tilt的三角函数型非线性方程组(6)，转化后的2次多项式型方程组的求解效率确实会提高。即使这样，转化后的求解方案也还是只能采用数值解法，还是不能满足在线定标的要求。

上述转化虽然没能完全解决问题，却为寻找未知参数ΔP和Tilt的闭式解带来了新的视角。将式(6)中含有sin(Tilt)和cos(Tilt)的项合并整理后，可变换为如下形式：

$\left(\begin{array}{cc}U{Y}_p\cos \left(\mathrm{\ΔP}\right)-U{X}_p\sin \left(\mathrm{\ΔP}\right)& Z_{p}U\\ V{X}_p\sin \left(\mathrm{\ΔP}\right)-V{Y}_p\cos \left(\mathrm{\ΔP}\right)-Z_p& Y_p\cos \left(\mathrm{\ΔP}\right)-V{Z}_p-X_p\sin \left(\mathrm{\ΔP}\right)\end{array}\right)\left(\begin{array}{c}\cos \left(\mathrm{Tilt}\right)\\ \sin \left(\mathrm{Tilt}\right)\end{array}\right)$

$=\left(\begin{array}{c}{X}_p\cos \left(\mathrm{\ΔP}\right)+Y_p\sin \left(\mathrm{\ΔP}\right)\\ 0\end{array}\right)---\left(7\right)$

式(7)可被看作是关于sin(Tilt)和cos(Tilt)的一个带参数ΔP的线性方程组。通过求解这个带参数ΔP的关于sin(Tilt)和cos(Tilt)的2阶线性方程组(7)，可以得到sin(Tilt)和cos(Tilt)的关于参数ΔP的闭式表达式。

再利用sin(Tilt)和cos(Tilt)要满足的平方和为1的约束，可以得到一个仅含有未知参数ΔP的方程。

若由这个仅含有未知参数ΔP的方程，可以得到未知参数ΔP的闭式解，则可将这个参数ΔP的闭式解代入sin(Tilt)和cos(Tilt)的关于参数ΔP的闭式表达式，就可以获得sin(Tilt)和cos(Tilt)的闭式解。从而得到Tilt的闭式解。

(2.2)求解带参数的线性方程组

在实际应用环境中，可以限制标识点不位于摄像机的正上方或正下方。因此，X_p＝Y_p＝0不会成立。

下面具体来讨论关于sin(Tilt)和cos(Tilt)的，带参数ΔP的线性方程组(7)的求解。

为了简化表达，记，

$A\left(\mathrm{\ΔP}\right)=\left(\begin{array}{cc}U{Y}_p\cos \left(\mathrm{\ΔP}\right)-U{X}_p\sin \left(\mathrm{\ΔP}\right)& Z_{p}U\\ V{X}_p\sin \left(\mathrm{\ΔP}\right)-V{Y}_p\cos \left(\mathrm{\ΔP}\right)-Z_p& Y_p\cos \left(\mathrm{\ΔP}\right)-V{Z}_p-X_p\sin \left(\mathrm{\ΔP}\right)\end{array}\right)$

$m\left(\mathrm{\ΔP}\right)=\left(\begin{array}{c}{X}_p\cos \left(\mathrm{\ΔP}\right)+Y_p\sin \left(\mathrm{\ΔP}\right)\\ 0\end{array}\right)$

则式(7)可转化为A(ΔP)W(Tilt)＝m(ΔP)；其中，A(ΔP)为关于sin(Tilt)和cos(Tilt)的线性方程组(7)中带参数ΔP的2阶系数矩阵，m(ΔP)为线性方程组(7)中带参数ΔP的2维向量，W(Tilt)＝(cos(Tilt)，sin(Tilt))^T是线性方程组(7)中待求解的2维变量向量。

记A(ΔP)的行列式为det(A(ΔP))，

det(A(ΔP))＝U[Z_p ²+(Y_pcos(ΔP)X_psin(ΔP))²]

下面讨论在各种情形下，方程组(7)的求解方法。

(2.2.1)det(A(ΔP))≠0

det(A(ΔP))≠0时，满足式(7)的sin(Tilt)和cos(Tilt)的闭式表达式为

$\cos \left(\mathrm{Tilt}\right)=\frac{(Y_p\cos \left(\mathrm{\ΔP}\right)-V{Z}_p-X_p\sin \left(\mathrm{\ΔP}\right))(X_p\cos \left(\mathrm{\ΔP}\right)+Y_p\sin \left(\mathrm{\ΔP}\right))}{U[{Z_p}^2+{(Y_p\cos \left(\mathrm{\ΔP}\right)-X_p\sin \left(\mathrm{\ΔP}\right))}^2]}$

$\sin \left(\mathrm{Tilt}\right)=\frac{(Z_p+V{Y}_p\cos \left(\mathrm{\ΔP}\right)-V{X}_p\sin \left(\mathrm{\ΔP}\right))(X_p\cos \left(\mathrm{\ΔP}\right)+Y_p\sin \left(\mathrm{\ΔP}\right))}{U[{Z_p}^2+{(Y_p\cos \left(\mathrm{\ΔP}\right)-X_p\sin \left(\mathrm{\Δp}\right))}^2]}---\left(8\right)$

将式(8)代入如下的三角函数等式约束sin²(Tilt)+cos²(Tilt)＝1，得到：

$[(V^2+1)Y_p^2-U^2(X_p^2+Z_p^2)]{\tan }^2\left(\mathrm{\ΔP}\right)+2X_p{Y}_p(U^2+V^2+1)\tan \left(\mathrm{\ΔP}\right)$

$+(V^2+1)X_p^2-U^2{Z}_p^2-U^2{Y}_p^2=0---\left(9\right);$

式(9)是一个关于tan(ΔP)的方程。

记 $a=(V^2+1)Y_p^2-U^2(X_p^2+Z_p^2),$ b＝2X_pY_p(U²+V²+1)， $c=(V^2+1)X_p^2-U^2{Z}_p^2-U^2{Y}_p^2$

式(9)可表示为，

atan²(ΔP)+btan(ΔP)+c＝0

a≠0时，式(9)是一个关于tan(ΔP)的二次方程。

记

$\mathrm{\Δ}=4U^2(X_p^2+Y_p^2+Z_p^2)[(X_p^2+Y_p^2)(V^2+1)-U^2{Z}_p^2]$

Δ≥0时，关于tan(ΔP)的二次方程(9)的二个闭式解可以表示为

$\tan \left(\mathrm{\ΔP}\right)=\frac{-X_p{Y}_p(U^2+V^2+1)+U\sqrt{(X_p^2+Y_p^2+Z_p^2)[(X_p^2+Y_p^2)(V^2+1)-U^2{Z}_p^2]}}{(V^2+1)Y_p^2-U^2(X_p^2+Z_p^2)}$

和

$\tan \left(\mathrm{\ΔP}\right)=\frac{-X_p{Y}_p(U^2+V^2+1)+U\sqrt{(X_p^2+Y_p^2+Z_p^2)[(X_p^2+Y_p^2)(V^2+1)-U^2{Z}_p^2]}}{(V^2+1)Y_p^2-U^2(X_p^2+Z_p^2)}$

这两个tan(ΔP)的闭式解中有一个是真解。

Tilt角的变化范围为[-78°，0°]，因此cos(Tilt)为正，sin(Tilt)为负，利用这个限制条件，可以从上面的tan(ΔP)的两个闭式解中选择真解。

具体来说，由于ΔP属于[-0.5°，0.5°]，因此，由tan(ΔP)的闭式解可得到sin(ΔP)和cos(ΔP)的如下形式的闭式解，

$\cos \left(\mathrm{\ΔP}\right)=\frac{1}{\sqrt{1+{\tan }^2\left(\mathrm{\ΔP}\right)}}$

$\sin \left(\mathrm{\ΔP}\right)=\frac{\tan \left(\mathrm{\ΔP}\right)}{\sqrt{1+{\tan }^2\left(\mathrm{\ΔP}\right)}}$

这样，根据tan(ΔP)的二个闭式解，可以确定出sin(ΔP)和cos(ΔP)的二组闭式解。将这二组闭式解代入到式(8)中，可以得到cos(Tilt)和sin(Tilt)的二组闭式解。其中，cos(Tilt)为正，sin(Tilt)为负的一组闭式解为真解。

对应于cos(Tilt)和sin(Tilt)的真解的那一个tan(ΔP)的闭式解为真解。

确定了tan(ΔP)的闭式解后，就能够确定cos(ΔP)和sin(ΔP)的闭式解。

这样，就得到了sin(ΔP)和cos(ΔP)，cos(Tilt)和sin(Tilt)的闭式解。从而得到ΔP和Tilt的闭式解。

综合上述求解过程，在

det(A(ΔP))＝U[Z_p ²+(Y_pcos(ΔP)-X_psin(ΔP))²]≠0，

$a=(V^2+1)Y_p^2-U^2(X_p^2+Z_p^2)\≠0$

$\mathrm{\Δ}=4U^2(X_p^2+Y_p^2+Z_p^2)[(X_p^2+Y_p^2)(V^2+1)-U^2{Z}_p^2]\≥0$

的条件下，先从式(7)得到含有未知参数sin(ΔP)和cos(ΔP)的，关于sin(Tilt)和cos(Tilt)的闭式解表达式(8)。

然后，利用三角函数等式，

sin²(Tilt)+cos²(Tilt)＝1

得到一个含有未知数tan(ΔP)的方程(9)。

由方程(9)解算出tan(ΔP)的闭式解，进而解算出sin(ΔP)和cos(ΔP)的闭式解。

最后，将sin(ΔP)和cos(ΔP)的闭式解代入到式(8)中，就可以解算出sin(Tilt)和cos(Tilt)的闭式解。

(2.2.2)det(A(ΔP))≠0并且a＝0

det(A(ΔP))＝U[Z_p ²+(Y_pcos(ΔP)X_psin(ΔP))²]≠0，

$a=(V^2+1)Y_p^2-U^2(X_p^2+Z_p^2)=0$

式(9)可表示为，

btan(ΔP)+c＝0

如果假定b＝2X_pY_p(U²+V²+1)≠0，则式(9)的解可以表示为

$\tan \left(\mathrm{\ΔP}\right)=\frac{U^2{Z}_p^2+U^2{Y}_p^2-(V^2+1)X_p^2}{2X_p{Y}_p(U^2+V^2+1)}---\left(10\right)$

实际上，若b＝2X_pY_p(U²+V²+1)＝0，则由式(9)，必有，

$c=(V^2+1)X_p^2-U^2{Z}_p^2-U^2{Y}_p^2=0$

从a＝0、b＝0和c＝0，可以推出X_p＝Y_p＝Z_p＝0。

根据前面的限定条件，这种情况不会出现，因此可以确定必有b≠0。

于是，det(A(ΔP))≠0并且a＝0时，式(10)总成立。

式(10)给出了tan(ΔP)的闭式解。

由tan(ΔP)的闭式解可以唯一的确定出sin(ΔP)和cos(ΔP)的闭式解。

将sin(ΔP)和cos(ΔP)的闭式解代入式(8)中，可以得到sin(Tilt)和cos(Tilt)的闭式解。

(2.2.3)det(A(ΔP))＝0

det(A(ΔP))＝U[Z_p ²+(Y_pcos(ΔP)-X_psin(ΔP))²]＝0，

若U不为0，则必有，

$\left\{\begin{array}{c}{Z}_p=0\\ Y_p\cos \left(\mathrm{\ΔP}\right)-X_p\sin \left(\mathrm{\ΔP}\right)=0\end{array}\right.,$

此时，式(7)可整理为，

Y_psin(ΔP)+X_pcos(ΔP)＝0。

这样，可以得到X_p＝Y_p＝Z_p＝0。

根据前面的限定条件，这种情况不会出现。

因此，det(A)＝0时，必有，

U＝0。

此时，式(7)可整理为，

$\left\{\begin{array}{c}{X}_p\cos \left(\mathrm{\ΔP}\right)+Y_p\sin \left(\mathrm{\ΔP}\right)=0\\ (V{X}_p\sin \left(\mathrm{\ΔP}\right)-V{Y}_p\cos \left(\mathrm{\ΔP}\right)-Z_p)\cos \left(\mathrm{Tilt}\right)+(Y_p\cos \left(\mathrm{\ΔP}\right)-X_p\sin \left(\mathrm{\ΔP}\right)-V{Z}_p)\sin \left(\mathrm{Tilt}\right)=0\end{array}\right.$

于是，tan(ΔP)的解为

$\tan \left(\mathrm{\ΔP}\right)=-\frac{X_p}{Y_p}---\left(11\right)$

tan(Tilt)的解为

$\tan \left(\mathrm{Tilt}\right)=\frac{Z_p+V{Y}_p\cos \left(\mathrm{\ΔP}\right)-V{X}_p\sin \left(\mathrm{\ΔP}\right)}{Y_p\cos \left(\mathrm{\ΔP}\right)-X_p\sin \left(\mathrm{\ΔP}\right)-V{Z}_p}---\left(12\right)$

式(11)给出了tan(ΔP)的闭式解，由tan(ΔP)的闭式解可以唯一的确定出如下形式的sin(ΔP)和cos(ΔP)的闭式解：

$\cos \left(\mathrm{\ΔP}\right)=\frac{1}{\sqrt{1+{\tan }^2\left(\mathrm{\ΔP}\right)}}$

$\sin \left(\mathrm{\ΔP}\right)=\frac{\tan \left(\mathrm{\ΔP}\right)}{\sqrt{1+{\tan }^2\left(\mathrm{\ΔP}\right)}}$

结合sin(ΔP)和cos(ΔP)的闭式解，式(12)给出了tan(Tilt)的闭式解。

由tan(Tilt)的闭式解，可以确定出如下形式的sin(Tilt)和cos(Tilt)的闭式解。

$\cos \left(\mathrm{Tilt}\right)=\frac{1}{\sqrt{1+{\tan }^2\left(\mathrm{Tilt}\right)}}$

$\sin \left(\mathrm{Tilt}\right)=\frac{\tan \left(\mathrm{Tilt}\right)}{\sqrt{1+{\tan }^2\left(\mathrm{Tilt}\right)}}$

综述所述，我们给出了各种情形下sin(ΔP)、cos(ΔP)、sin(Tilt)和cos(Tilt)的闭式解。

由sin(ΔP)、cos(ΔP)、sin(Tilt)和cos(Tilt)的闭式解，可以得到ΔP和Tilt的闭式解。

由ΔP和P₀，可以得到Pan的闭式解。

由Pan和Tilt的闭式解，可以得到式(1)～(3)中的旋转变换矩阵R的闭式解，完成PTZ摄像机的单点在线标定。

如前所述，在户外恶劣的陌生环境中，有时无法获得足够数目的标识点的三维物理坐标信息。极端情形下，只有一个标识点的三维物理坐标信息可被用于Pan和Tilt的标定。在这种仅有一个标识点且需要在线标定的应用环境中，现有的标定方法都不适用。我们提出的这种PTZ摄像机单点在线标定方法是唯一选择。

这种标定方法在求解过程中，巧妙的将复杂费时的非线性方程组的求解问题转化为带参数的线性方程组的求解问题，并且在各种情形下给出了sin(ΔP)、cos(ΔP)、sin(Tilt)和cos(Tilt)的闭式解。

这种单点标定方法不仅是在只有一个标识点，且需要在线标定的应用环境中的唯一选择。而且，实验表明，采用这种单点标定方法得到的PTZ摄像机的水平转角Pan和垂直转角Tilt的绝对误差也非常低。

为了使用PTZ摄像机的单点在线标定方法，需要给定摄像机的内参数矩阵K、摄像机在基准坐标系中的三维物理坐标(x_c，y_c，z_c)^T、一个标识点在基准坐标系中的三维物理坐标(x，y，z)^T、该标识点在图像上所对应像点的像素坐标(u，v)^T、PTZ摄像机输出的Pan角度的整数形式的粗精度初始参考值P₀，Tilt需属于[-78°，0°]、该标识点需要不位于摄像机的正上方或正下方。满足这些条件后，这种PTZ摄像机的单点在线标定方法就能够在线提供PTZ摄像机的水平转角Pan和垂直转角Tilt的闭式解。

如图3所示，PTZ摄像机的单点标定方法的步骤具体如下：

S1：根据标识点在基准坐标系中的三维物理坐标和与所述标识点对应的像点在图像上的二维像素坐标获得标定的关系模型；

S2：采用线性方法计算所述标定的关系模型并获得水平转角和垂直转角。

其中步骤S2的子流程如图4所示，具体包括：

S21：将标定的关系模型转化成带水平转角参数并关于垂直转角的正、余弦的线性关系模型；

S22：根据所述线性的关系模型获得垂直转角的正、余弦的闭式表达式；

S23：根据垂直转角的正、余弦闭式表达式并利用三角函数约束得到水平转角的正切的闭式表达式；

S24：根据所述水平转角的正切闭式表达式和所述垂直转角的正、余弦闭式表达式得到水平转角和垂直转角。

为了更进一步说明本发明实施例提供的基于二自由度摄像机的单点标定方法，现通过仿真实验对单点标定的准确度做一个定量的描述。首先我们给定PTZ摄像机的内参数矩阵K，

$K=\left(\begin{array}{ccc}560.0& 0& 512.0\\ 0& 560.0& 512.0\\ 0& 0& 1\end{array}\right)$

在以下实验中，PTZ摄像机的内参数矩阵K固定且不变；设定PTZ摄像机在基准坐标系中的三维物理坐标(x_c，y_c，z_c)^T；设定标识点在基准坐标系中的三维物理坐标(x，y，z)^T。

如表1所示，设定10组不同的PTZ摄像机的Pan角和Tilt角真值。将Pan角和Tilt角代入式(2)，得到旋转矩阵R的真值。

表1设定的10组的PTZ摄像机的Pan角和Tilt角

编号	Pan(度)	Tilt(度)	编号	Pan(度)	Tilt(度)
						1	0.3	-0.3	6	180.56	-35.27
2	30.72	-7.82	7	228.68	-42.89
						3	79.63	-14.29	8	270.61	-50.21
4	90.35	-21.85	9	332.57	-62.19
						5	142.28	-28.45	10	359.81	-77.82

将矩阵K、坐标(x_c，y_c，z_c)^T、坐标(x，y，z)^T和矩阵R代入式(1)中，得到标识点所对应像点在图像中的像素坐标的精确值。

由于在实际应用中，获取的标识点所对应像点的像素坐标可能存在1个像素的误差。因此，我们分别设定标识点所对应像点的像素坐标为：(1)像素坐标的精确值；(2)取整后的像素坐标；(3)取整后像素坐标的四邻域像素坐标。并在以上标识点所对应像点的像素坐标(u，v)^T的不同设定方式下，对单点标定方法进行测试。这样，依据设定的矩阵K、坐标(x_c，y_c，z_c)^T、坐标(x，y，z)^T、坐标(u，v)^T，采用单点标定算法，我们可以解算出水平转角Pan和垂直转角Tilt。

将解算出的水平转角Pan和垂直转角Tilt与设定的水平转角Pan和垂直转角Tilt的真值进行比对，计算出单点标定算法解算的水平转角Pan和垂直转角Tilt的绝对误差，比对结果如图5和图6所示。图5和图6分别绘制了在设定标识点所对应像点的像素坐标为：精确像素坐标，取整后像素坐标和取整后像素坐标的四邻域像素坐标时，由单点标定方法解算的Pan和Tilt的绝对误差，以及不采用单点定标方法直接采用整数形式的Pan和Tilt的输出值时的绝对误差。从图5和图6可以看出，采用精确像素坐标，取整后像素坐标和取整后像素坐标的四邻域像素坐标由单点标定方法解算的Pan和Tilt的绝对误差要比不采用单点标定直接采用整数形式的Pan和Tilt的输出值的绝对误差要低。采用单点标定算法提高了摄像机标定的准确度。如果设定标识点所对应像点的像素坐标为精确像素坐标，采用单定定标方法解算出的Pan和Tilt绝对误差几乎为零。设定标识点所对应像点的像素坐标为取整后的像素坐标，采用单定定标方法解算出的Pan和Tilt绝对误差也比较小。设定标识点所对应像点的像素坐标为取整后的四邻域像素坐标，采用单定定标方法解算出的Pan和Tilt绝对误差虽然有所增加。但是比不采用单点定标方法而直接采用整数形式的Pan和Tilt的输出值的绝对误差要低。也就是说，即便给定的标识点的像素坐标可能存在误差，也要比直接采用整数形式的Pan和Tilt的输出值而不标定的效果要好。

通过以上的实验分析可以看出，这种PTZ摄像机的单点在线标定方法不仅是在仅有一个标识点的应用环境中的唯一选择，而且性能也很好。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于二自由度摄像机的单点标定方法，其特征在于，包括下述步骤：

S1：根据标识点在基准坐标系中的三维物理坐标和与所述标识点对应的像点在图像上的二维像素坐标获得标定的关系模型；

S2：采用线性方法计算所述标定的关系模型并获得水平转角和垂直转角。

2.如权利要求1所述的单点标定方法，其特征在于，所述标定的关系模型具体为： $Z\left(\begin{array}{c}U\\ 1\\ -V\end{array}\right)=\left(\begin{array}{ccc}\cos \left(\mathrm{\ΔP}\right)& \sin \left(\mathrm{\ΔP}\right)& 0\\ -\sin \left(\mathrm{\ΔP}\right)\cos \left(\mathrm{Tilt}\right)& \cos \left(\mathrm{\ΔP}\right)\cos \left(\mathrm{Tilt}\right)& \sin \left(\mathrm{Tilt}\right)\\ \sin \left(\mathrm{\ΔP}\right)\sin \left(\mathrm{Tilt}\right)& -\cos \left(\mathrm{\ΔP}\right)\sin \left(\mathrm{Tilt}\right)& \cos \left(\mathrm{Tilt}\right)\end{array}\right)\left(\begin{array}{c}{X}_p\\ Y_p\\ Z_p\end{array}\right);$

U表示标识点在焦距为1时的正则化像平面上所对应像点的像素列坐标，V表示标识点在焦距为1时的正则化像平面上所对应像点的像素行坐标，Z表示标识点P在像空坐标系C-XYZ中的第三维坐标，(X_p，Y_p，Z_p)^T表示标识点在临时坐标系C-x’y’z中的三维物理坐标；ΔP表示水平转角Pan与初始参考值P₀之间的差值，Tilt表示垂直转角。

3.如权利要求2所述的单点标定方法，其特征在于，所述临时坐标系C-x’y’z的获取方法包括：

将基准坐标系O-xyz平移得到过渡临时坐标系C-xyz；

再将过渡临时坐标系C-xyz绕其Cz轴正向右手螺旋旋转角度P₀后得到所述临时坐标系C-x’y’z。

4.如权利要求1或2所述的单点标定方法，其特征在于，所述步骤S2具体包括：

S21：将所述标定的关系模型转化成带水平转角参数并关于垂直转角的正、余弦线性关系模型；

S22：根据所述线性的关系模型获得垂直转角的正、余弦闭式表达式；

S23：根据垂直转角的正、余弦闭式表达式并利用三角函数约束得到水平转角的正切闭式表达式；

S24：根据所述水平转角的闭式表达式和所述垂直转角的闭式表达式得到水平转角和垂直转角。

5.如权利要求4所述的单点标定方法，其特征在于，所述线性关系模型具体为：

$\left(\begin{array}{cc}U{Y}_p\cos \left(\mathrm{\ΔP}\right)-U{X}_p\sin \left(\mathrm{\ΔP}\right)& Z_{p}U\\ V{X}_p\sin \left(\mathrm{\ΔP}\right)-V{Y}_p\cos \left(\mathrm{\ΔP}\right)-Z_p& Y_p\cos \left(\mathrm{\ΔP}\right)-V{Z}_p-X_p\sin \left(\mathrm{\ΔP}\right)\end{array}\right)\left(\begin{array}{c}\cos \left(\mathrm{Tilt}\right)\\ \sin \left(\mathrm{Tilt}\right)\end{array}\right)$

$=\left(\begin{array}{c}{X}_p\cos \left(\mathrm{\ΔP}\right)+Y_p\sin \left(\mathrm{\ΔP}\right)\\ 0\end{array}\right).$

6.如权利要求4或5所述的单点标定方法，其特征在于，当线性关系模型的系数矩阵行列式不等于零时，垂直转角的正、余弦闭式表达式为：

$\cos \left(\mathrm{Tilt}\right)=\frac{(Y_p\cos \left(\mathrm{\ΔP}\right)-V{Z}_p-X_p\sin \left(\mathrm{\ΔP}\right))(X_p\cos \left(\mathrm{\ΔP}\right)+Y_p\sin \left(\mathrm{\ΔP}\right))}{U[{Z_p}^2+{(Y_p\cos \left(\mathrm{\ΔP}\right)-X_p\sin \left(\mathrm{\ΔP}\right))}^2]}$ 和

$\sin \left(\mathrm{Tilt}\right)=\frac{(Z_p+V{Y}_p\cos \left(\mathrm{\ΔP}\right)-V{X}_p\sin \left(\mathrm{\ΔP}\right))(X_p\cos \left(\mathrm{\ΔP}\right)+Y_p\sin \left(\mathrm{\ΔP}\right))}{U[{Z_p}^2+{(Y_p\cos \left(\mathrm{\ΔP}\right)-X_p\sin \left(\mathrm{\Δp}\right))}^2]};$

水平转角的正切闭式表达式为：

$\tan \left(\mathrm{\ΔP}\right)=\frac{-X_p{Y}_p(U^2+V^2+1)\text{\±U}\sqrt{(X_p^2+Y_p^2+Z_p^2)[(X_p^2+Y_p^2)(V^2+1)-U^2{Z}_p^2]}}{(V^2+1)Y_p^2-U^2(X_p^2+Z_p^2)}$ 或

$\tan \left(\mathrm{\ΔP}\right)=\frac{U^2{Z}_p^2+U^2{Y}_p^2-(V^2+1)X_p^2}{2X_p{Y}_p(U^2+V^2+1)}.$

7.如权利要求4或5所述的单点标定方法，其特征在于，当线性关系模型的系数矩阵行列式等于零时，所述垂直转角的正切闭式表达式为 $\tan \left(\mathrm{Tilt}\right)=\frac{Z_p+V{Y}_p\cos \left(\mathrm{\ΔP}\right)-V{X}_p\sin \left(\mathrm{\ΔP}\right)}{Y_p\cos \left(\mathrm{\ΔP}\right)-X_p\sin \left(\mathrm{\ΔP}\right)-V{Z}_p};$ 所述水平转角的正切闭式表达式为

$\tan \left(\mathrm{\ΔP}\right)=-\frac{X_p}{Y_p}.$

Images