CN103454276B - 一种基于动态序列图像的织物形态风格评价方法 - Google Patents

Abstract

本发明涉及一种织物形态风格评价方法，属于纺织领域。一种基于动态序列图像的织物形态风格评价方法，其特征在于利用CCD数码相机对织物面料的运动过程进行动态拍摄，得到织物运动的多帧动态序列视图；从同一帧的所有拍摄的视图上提取特征点；将对应视图上的特征点进行立体匹配并计算匹配点的三维坐标；从获得的特征点的三维坐标计算织物的基本运动物理量；根据织物的基本运动物理量及运动特征提取织物动态特征参数；对织物的物理性能进行测量；对这些测量的物理量进行因子分析；计算动态特征参数与织物物理性能间的关联关系，结合基于物理性能测量中提取形态风格评价因子，得到基于动态序列图像的形态风格评价因子，完成织物形态的评价。

Description

一种基于动态序列图像的织物形态风格评价方法

技术领域

本发明涉及一种织物形态风格评价方法，尤其涉及一种基于动态序列图像的织物形态风格评价方法。

背景技术

织物的形态特征是指织物在赋形中形成的线条形态刺激视觉器官后产生的一种涉及心理映射的感觉效果。它可以通过织物的刚柔性、赋形性和形态稳定性等物理概念间接加以表达，即服装织物的物理力学性能（如拉伸、剪切、弯曲等）影响着服装的外观。

70年代,以京都大学高分子化学系教授川端季雄为首的织物风格评价与标准化委员会（HESC）开发了川端风格仪（KES-F）系统，该系统基于尽可能全面地反映织物特性的出发点，从织物的力学性能如弯曲刚度、剪切刚度等以及其他一些复合性能（如由表面摩擦和几何学的粗糙度与力学性能混合感觉到的表面光滑性），来评价织物的风格感觉特性。他们通过测试拉伸、压缩、剪切、弯曲和表面性能五项基本力学性能中的16项物理指标，再加上单位面积质量来推断织物的手感风格。该系统在大量工作的基础上，将不同用途织物的风格分解成若干个基本风格，并将综合风格和基本风格量化，分别建立织物物理力学量和基本风格值之间、基本风格值和综合风格值之间的回归方程式，并根据织物的基本风格来决定服装的用途。

Postle等人开发出了可用于织物外观、手感和机械性能的简便测试与快速评价的FAST系统。与KES-F系统不同的是，FAST系统基于简单测量的出发点，有选择性地测量织物在小应力、小变形条件下的压缩、弯曲、拉伸、剪切四项基本力学性能和尺寸稳定性，其力学性能实质上只是简单地测试了织物在一些特征性的小应力条件下的变形，而且未像KES-F那样测试织物的变形回复过程。

已有的织物风格评价的研究现状表明，现有的评价标准不能真正从客观上全面、方便地评价织物的形态风格特征，主要表现在以下两个方面：第一，在织物形态风格的测定和评价中，三维动态风格参数的分析、处理较少，各评价指标中的测量数据多由织物静态时测得；第二，分指标多次测量再综合评价的方法增加了评价工作的繁复性，测量时得出的数据受主观因素的影响较大，且有些指标对实际生产检测的意义不大。

我国作为一个纺织大国，近年来纺织行业取得了飞速的发展，在服装企业中，织物形态风格的优劣对于织物的选用和开发、服装设计具有至关重要的作用。如果能有一种完善可行的织物形态风格的评价标准，就可以预测或区分各织物和服装的动态外观效果，从而可大大提高织物和服装产品在设计、开发、销售等各环节的效率。

目前用仪器客观地测试织物的物理机械特性已不是难事，但如何把仪器测试值与人的感官评定结果联系起来，即寻找织物物理量与风格心理量间的关系表达式，一直是研究的重点。尽管研究取得了很大进展，但一直存在一个最基本的问题，那就是：在以往的评价中，所采取的测量方法主观性比较大，获得的测量数据受人的主观因素比较大。因此，在织物的形态风格评价领域，对织物的形态采取怎样的检测方法，进而选取表征织物形态的特征参数是关键。因此亟待发展一种非接触式的测量方法去科学评价织物的风格。多目立体视觉是计算机立体视觉中的重要分支，同时是近年来计算机视觉中研究的热点。基于多目立体视觉的图像配准技术是根据人在大脑中将左右眼图像融合成立体图像的原理，使用标定后的若干个摄像头进行同步拍摄获得同一物体在相同条件下不同位置的图像，然后根据所获取的二维图像的视差和像素点灰度信息进行计算，获取像素点的第三维深度信息。关于双目立体视觉的应用层出不穷，如RyanWhiteRyanWhite采用多台高速摄像机对印有特殊图案织物进行运动捕捉，用近似推理的思想挖掘局部领域信息与整体应力信息，以此来确定褶皱在织物上堆积的位置，运用该方法，实现了手臂弯曲时在肘部形成的织物褶皱形状、裤子快速移动时的形态、服装突然被触及时的相应织物部位形状的三维重建；DanielVlasic、Vogiatzis、Hernandez也运用多角度的广度立体视差法对抓取服装的运动状态进行了研究，利用特征点位置在织物三维重建的表面建立三角网格，再利用插值对表面进行光滑处理，较好地获得了织物在虚拟世界中的展示效果。那么，利用多角度的高速摄像机可以从整体上客观、精确地把握织物的运动情况，为科学地描述织物动态外观形态提供了参考。本发明正是围绕双目立体视觉对织物的动态捕捉而展开的。

发明内容

本发明所要解决的技术问题是提供一种基于动态序列图像的织物形态风格评价方法，解决对于织物的形态风格评价中，采取的测量方法和检测评价方法主观性过大，而且需要分指标进行多次测量再综合评价，工作繁复且数据不精确的问题。

技术方案

一种基于动态序列图像的织物形态风格评价方法，其特征在于：包括利用设置的多台CCD数码相机对织物面料的运动过程进行动态拍摄，得到织物运动的多帧动态序列视图的步骤；利用计算机视觉处理技术从同一帧的所有拍摄的视图上提取特征点的步骤；将对应视图上的特征点进行立体匹配并计算匹配点在世界坐标系中的三维坐标的步骤；从获得的特征点的三维坐标计算织物的基本运动物理量的步骤，所述基本运动物理量包括各特征点在X,Y,Z方向上的位移、速度、加速度；根据织物的所述基本运动物理量及运动特征提取织物动态特征参数的步骤，所述织物动态特征参数包括织物面料在Z方向上能达到的最大累积位移、运动过程中织物表现出来的“波”形曲线在最凹处时的曲率以及“波”形曲线在运动过程中的保形率；采用传统的测量方法对织物的物理性能进行测量的步骤，所述的物理性能包括织物结构参数与物理机械性能；然后，对这些测量的物理性能进行因子分析的步骤,提取计算与织物形态风格相关的3个主因子，所述3个主因子分别为悬垂投影因子、悬垂美感因子、波均匀度因子；计算所述动态特征参数与织物物理性能间的关联关系，结合对测量的物理性能进行因子分析的步骤中提取计算的形态风格评价因子，得到基于动态序列图像的3个形态风格评价主因子，完成织物形态在视觉方面的形态风格评价。

进一步，所述传统的测量方法包括使用FAST织物风格仪、静态悬垂仪和密度镜对对织物的物理性能进行测量。

进一步，所述提取特征点的步骤中，采用Canny边缘检测算法和Harris角点检测算法结合进行特征点的提取，再采用SIFT特征点匹配算法和双目视觉视差原理获得各特征点在世界坐标系中的三维坐标。

进一步，所述得到的基于动态序列图像的3个形态风格评价主因子，完成织物形态在视觉方面的评价，其中所述3个主因子的计算式如下：

进一步，将待评测的多种不同的织物面料进行裁剪，裁剪后将各织物面料的运动过程进行单个的多帧动态视图拍摄，提取出织物动态特征参数，依据参数进行对比和评价，所述待评测的多种不同的织物面料采用同一条件进行拍摄和织物动态特征提取。

所述同一条件包括各织物面料平面尺寸的大小相同、触发面料运动的机制条件相同。

有益效果

本发明的方法分别基于织物静态物理性能和动态序列图像两方面对织物的风格进行了评价，只要对织物进行一次物理性能测量，然后从动态序列图像中分析计算得到各运动特征参数，建立织物的动态特征参数与织物静态物理性能间的关联关系，根据关联关系得到基于动态序列图像计算的动态特征参数的三个体现织物形态风格的主因子，即可完成对织物的风格评价。本发明解决了人眼视觉（即风格心理量）与静态物理量的联通关系，提出了三个无论在静态还是动态都能具体表征织物形态的特征参数（即三个形态风格评价主因子），实现了快速、简便、客观地完成对织物的风格评价。

附图说明

附图1为两台相对位置任意的（镜头均朝向被拍摄目标）高速摄像机拍摄目标P的示意图；

附图2为摄像机透视投影模型；

附图3为图像坐标系；

附图4为左相机部分特征点的图像坐标；

附图5为右相机部分特征点的图像坐标；

附图6为各特征点的三维世界坐标；

附图7为各特征点在左相机中的计算可靠性；

附图8为各特征点在右相机中的计算可靠性；

附图9为静态悬垂仪测出的投影图，其中1为波峰半径R_peak；2为波谷半径R_valley；3为波峰间距D_peaks；4为波谷间距D_valleys；

附图10为涤纶仿丝软缎、真丝软缎、真丝绸、顺纡绉、平纹全棉、涤塔夫这六种织物的一个运动周期中的部分视图；

附图11～16分别为这涤纶仿丝软缎、真丝软缎、真丝绸、顺纡绉、平纹全棉、涤塔夫这六种织物的动态特征参数表；

附图17为涤纶仿丝软缎、真丝软缎、真丝绸、顺纡绉、平纹全棉、涤塔夫这六种织物的结构参数及物理机械性能；

附图18为因子分析中的解释的总方差；

附图19为各因子的最终得分矩阵；

附图20为各织物物理性能与Z方向上的最大累积位移的数量关系；

附图21为各织物物理性能与曲率ρ的数量关系；

附图22为各织物物理性能与曲线保形率α的数量关系。

具体实施方式

下面结合具体实施例和附图，进一步阐述本发明。

为了对织物的形态风格评价更客观，本发明提出一种基于动态序列图像的织物形态特征提取和织物形态风格评价方法，同时提出一些织物的动态特征参数和织物形态风格评价参数，从而依据这些参数进行客观比较和评价。

以附图1所示的两台单目高速摄像机作为实施例，即左右各设置一台来解析织物形态特征提取和织物形态风格评价方法，具体的步骤为：①利用设置的多台CCD数码相机对织物的运动过程进行动态拍摄，得到织物运动的多帧动态序列视图的步骤；②利用计算机视觉处理技术从同一帧的所有拍摄的视图上提取特征点的步骤；③将对应视图上的特征点进行立体匹配并计算匹配点在世界坐标系中的三维坐标的步骤；④从获得的特征点的三维坐标计算织物的基本运动物理量的步骤，所述基本运动物理量包括各特征点在X,Y,Z方向上的位移、速度、加速度；⑤根据织物的所述基本运动物理量及运动特征提取织物动态特征参数的步骤，所述织物动态特征参数包括织物在Z方向上能达到的最大累积位移、运动过程中织物表现出来的“波”形曲线在最凹处时的曲率以及“波”形曲线在运动过程中的保形率；⑥采用传统的测量方法对织物的物理性能进行测量，所述的物理性能包括织物结构参数与物理机械性能；⑦对这些测量的物理量进行因子分析,提取与织物形态风格相关的3个主因子，所述3个主因子分别为悬垂投影因子、悬垂美感因子、波均匀度因子；⑧利用SPSS数据分析软件计算动态特征参数与织物物理性能间的关联关系，结合在基于织物物理性能中提取的形态风格评价因子，得到基于动态序列图像的形态风格评价因子，从而完成织物形态在视觉方面的评价。

在第①步骤获得织物的动态序列图像之前，需要对两台高速摄像机分别进行标定并标定两摄像机的内外参数的关系，标定方法如下：

对单个摄像机的标定过程就是确定摄像机的几何与光学参数，包括M_internal（包括x方向尺度因子y方向尺度因子摄像机光轴与图像平面的交点(u₀,v₀)，畸变因子η）、（旋转矩阵R的3个偏转角α，β，γ和平移向量T的3个参数），进而获得摄像机中心点相对于世界坐标系的位置。摄像机标定研究的基本模型为针孔成像模型，然而，在很多情况下，这种线性模型不能精确地描述摄像机成像的几何关系。因此，在考虑线性或非线性的畸变偏差后，利用校正后的模型进行三维重建才能得到更为精确的位置信息。在普遍情况下，只考虑径向畸变系数，并假设k＝k₁＝k₂,其数学模型为：

φ_rx＝x[k₁(x²+y²)+k₂(x²+y²)²]

φ_rx＝y[k₁(x²+y²)+k₂(x²+y²)²]

其中，k₁，k₂为径向畸变系数。

经典的摄像机标定方法是利用已知物体上（如棋盘格标定板）的某些点的三维坐标和对应的图像点坐标，计算摄像机的内外部参数。经典的标定方法主要有透视变换法、张正友法、直接线性变换法、两步法和双平面法。在本发明中，采用了经典的张正友标定方法。

i单应性矩阵的计算

在标定过程中，为简化计算，设定模板（黑白棋盘格）所在平面为世界坐标系的Z_w＝0平面。用r_i表示R的第i列矢量，那么对于模板上的任一点都有

$z_p\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{ccc}\frac{1}{\mathrm{dX}}& 0& u_0\\ 0& \frac{1}{\mathrm{dY}}& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}f& 0& 0& 0\\ 0& f& 0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cc}R& T\\ 0& 1\end{array}\right]\left[\begin{array}{c}{X}_w\\ Y_w\\ Z_w\\ 1\end{array}\right]$

$=F\left[\begin{array}{ccc}{r}_1& r_2& T\end{array}\right]\left[\begin{array}{c}{X}_w\\ Y_w\\ 1\end{array}\right]{\mathrm{MW}}_{\mathrm{XYZ}}$

则模板平面上的点和它的像点之间的映射可由矩阵确定，矩阵M成为单应性矩阵。设 $M=\left[\begin{array}{ccc}{m}_{11}& m_{12}& m_{13}\\ m_{21}& m_{22}& m_{23}\\ m_{31}& m_{32}& m_{33}\end{array}\right],$ 可得 $\left\{\begin{array}{c}{z}_{p}u=m_{11}{X}_w+m_{12}{Y}_w+m_{13}\\ z_{p}v=m_{21}{X}_w+m_{22}{Y}_w+m_{23}\\ z_p=m_{31}{X}_w+m_{32}{Y}_w+1\end{array}\right.$ 整理得

$\left[\begin{array}{cccccccc}{X}_w& Y_w& 1& 0& 0& 0& {-\mathrm{uX}}_w& {-\mathrm{uY}}_w\\ 0& 0& 0& X_w& Y_w& 1& {-\mathrm{vX}}_w& {-\mathrm{vY}}_w\end{array}\right]m=\left[\begin{array}{c}u\\ v\end{array}\right]$

其中m＝[m₁₁m₁₂m₁₃m₂₁m₂₂m₂₃m₃₁m₃₂m₃₃]^T

令M＝[m₁m₂m₃]，则

[m₁m₂m₃]＝λF[r₁r₂T]

因为r₁和r₂是单位正交矢量，即r₁ ^Tr₁＝r₂ ^Tr₂＝１，且r₁ ^Tr₂＝0，所以

m₁ ^TF^-TF^-1m₂＝０

m₁ ^TF^-TF^-1m₁＝m₂ ^TF^-TF^-1m₂

上两式为摄像机内部参数求解的两个约束条件式。

令

$B=F^{-T}{F}^{-1}=\left[\begin{array}{ccc}{b}_{11}& b_{12}& b_{13}\\ b_{21}& b_{22}& b_{23}\\ b_{31}& b_{32}& b_{33}\end{array}\right]=$

$\left[\begin{array}{ccc}\frac{1}{{f_u}^2}& -\frac{\mathrm{\η}}{{f_u}^2{f}_v}& \frac{v_0\mathrm{\η}-u_0{f}_v}{{f_u}^2{f}_v}\\ -\frac{\mathrm{\η}}{{f_u}^2{f}_v}& \frac{{\mathrm{\η}}^2}{{f_u}^2{f}_v}+\frac{1}{{f_v}^2}& -\frac{\mathrm{\η}(v_0\mathrm{\η}-u_0{f}_v)}{{f_u}^2{f}_v}-\frac{v_0}{{f_v}^2}\\ \frac{v_0\mathrm{\η}-u_0{f}_v}{{f_u}^2{f}_v}& -\frac{\mathrm{\η}(v_0\mathrm{\η}-u_0{f}_v)}{{f_u}^2{f}_v}-\frac{v_0}{{f_v}^2}& -\frac{\mathrm{\η}(v_0\mathrm{\η}-u_0{f}_v)}{{f_u}^2{f}_v}+\frac{v_0}{{f_v}^2}+1\end{array}\right]$

从上式可得，B为对称矩阵，定义6维矢量b＝[b₁₁b₁₂b₁₃b₂₁b₂₂b₂₃]，则

有m_i ^TBm_j＝v_ij ^Tb

其中，

m_i＝[m_i1m_i2m_i2]^T，

v_ij＝[m_i1m_j1m_i1m_j2+m_i2m_j1m_i2m_j2m_i3m_j1+m_i1m_j3m_i3m_j2+m_i2m_j3m_i3m_j3]

根据两个约束条件式，可得关于矢量b的两个齐次方程，即

$\left[\begin{array}{c}{v_{12}}^T\\ {(v_{11}-v_{22})}^T\end{array}\right]b=0$

若有n幅平面模板图像，根据上式，可得到如下线性方程组为：Vb＝0求得b即可获得矩阵B，进而求得摄像机的各内部参数，即

$\left\{\begin{array}{c}k=b_{33}-[{b_{13}}^2+v_0(b_{12}{b}_{13}-b_{11}{b}_{23})]/b_{11}\\ v_0=(b_{12}{b}_{12}-b_{11}{b}_{23})/(b_{11}{b}_{22}-{b_{12}}^2)\\ u_0={\mathrm{\ηv}}_0/f_u-b_{13}{f_u}^2/k\\ f_v=\sqrt{{\mathrm{kb}}_{11}(b_{11}{b}_{22}-{b_{12}}^2)}\\ f_u=\sqrt{k/b_{11}}\\ \mathrm{\η}=-b_{12}{f_u}^2{f}_v/k\end{array}\right.$

利用F与单应性矩阵M对每个不同姿态的摄像机外参数可按如下公式求得：

r₁＝λF^-1m₁，r₂＝λF^-1m₂，r₃＝r₁r₂，T＝λF^-1m₃,λ＝1/||F^-1m₁||＝1/||F^-1m₂||。由于噪声的存在，R＝[r₁r₂r₃]一般不满足旋转矩阵的性质，因此，应采用最大似然估计对R进行精确求解。

ii最大似然估计

设有n幅关于模板平面的图像，每幅图像中有m个标定点，每个标定点的数据都被独立同分布的噪声污染。建立如下的目标函数：

$\min \underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\left|\right|p_{\mathrm{ij}}-\hat{p}(F,R_i,T_i,M_j)\left|\right|}^2$

式中，p_ij为第i幅图像中第j个标定点的像点坐标；R_i、T_i分别为第i幅图像坐标系的旋转矩阵和平移矩阵，是第i幅图像中的第j标定点的世界坐标的M_j求单应性矩阵而求得的像点坐标。

iii径向畸变系数的计算

在标定过程中，必须考虑摄像机镜头畸变，实际上，当摄像机镜头远离图像中心较远处时，会产生较大的几何畸变。在此过程中，主要考虑一次与二次径向畸变，其他畸变因素忽略。摄像机镜头径向畸变模型为

$\left\{\begin{array}{c}{x_p}^{\′}={x_p}^{*}+{x_p}^{*}[k_1({x_p}^{*2}+{y_p}^{*2})+k_2{({x_p}^{*2}+{y_p}^{*2})}^2]\\ {y_p}^{\′}={y_p}^{*}+{y_p}^{*}[k_1({x_p}^{*2}+{y_p}^{*2})+k_2{({x_p}^{*2}+{y_p}^{*2})}^2]\end{array}\right.$

根据像平面和图像坐标的转换关系 $\left\{\begin{array}{c}u=\frac{{x_p}^{\′}}{\mathrm{dX}}+{{\mathrm{\ηy}}_p}^{\′}+u_0\\ v=\frac{{y_p}^{\′}}{\mathrm{dY}}+v_0\end{array}\right.,$ 可得

$\left\{\begin{array}{c}u=u^{\′}+(u^{\′}+u_0)[k_1({(u^{\′}+u_0)}^2+{(v^{\′}+v_0)}^2)+k_2{({(u^{\′}+u_0)}^2+{(v^{\′}+v_0)}^2)}^2]\\ v=v^{\′}+(v^{\′}+v_0)[k_1({(u^{\′}+u_0)}^2+{(v^{\′}+v_0)}^2)+k_2{({(u^{\′}+u_0)}^2+{(v^{\′}+v_0)}^2)}^2]\end{array}\right.$

(u,v)为无畸变影响的理想的图像像素点坐标，(u',v')为图像实际像素点坐标。那么，对于给定的n幅标定图像，每幅m个点，共mn个标定点，可以得到2mn个方程，即可以求得k₁与k₂，就可以求得相机内外参数的最优解。

iv.双摄像机立体系统标定

首先利用单摄像机标定方法对双目立体视觉中左右摄像机的外部参数R₁、T₁和R₂、T₂，则双视觉系统两个摄像机之间的相对几何位置关系可用下式表示：

$\left\{\begin{array}{c}R=R_2{R_1}^{-1}\\ T=T_2-R_2{R_1}^{-1}{T}_1\end{array}\right.$

在第②步骤提取特征点中，在matlab中，先采用Canny算子进行边缘检测，其检测步骤为：I.用标准偏差γ的高斯滤波器对图像进行平滑处理，用以抑制噪声；II.对每一像素点，计算局部梯度和边缘方向α(x,y)＝arctan(G_y/G_x)。梯度方向上强度最大的点即为边缘点；III.在上一步中检测到的边缘点会导致梯度的幅度图像中出现脊。对脊像素用两个阈值T₁和T₂做阈值处理，其中T₁＜T₂。值大于T₂的脊像素为强边缘像素，值处于T₁和T₂之间的脊像素为弱边缘像素；IV.连接边缘。

再采用Harris算子提取物体的角点信息，Harris检测算法的数学表达式为： $E(u,v)=\underset{x,y}{\mathrm{\Σ}}w(x,y){[I(x+u,y+v)-I(x,y)]}^2,$ 将图像窗口平移[u,v]产生灰度变化E(u,v)。其中w(x,y)为窗口函数，I(x,y)为图像灰度，I(x+u,y+v)为平移后的图像灰度。其检测算法的步骤为：I.计算图像像素在水平和垂直方向上的梯度，以及两者的乘积，得到如下矩阵： $M=\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& {I_y}^2\end{array}\right];$ II.对图像进行高斯滤波，得到新的M。离散二维零均值高斯函数为III.计算原图像上对应的每个像素点的兴趣值R：IV.选取局部范围内的极大兴趣值对应的像素点即为特征点；V.结合边缘检测获取的边缘点，设定阈值，选取一定量的角点。

在第③步骤计算特征点的世界坐标中，需要解决摄像机透视投影原理。如附图2所示，摄像机透视投影模型，即摄像机成像模型，是通过成像透镜将三维场景信息投射到摄像机二维平面上，且该投影可以用成像变换描述。

其中点O是摄像机的光心，摄像机坐标系的x轴与y轴与图像坐标系的X轴与Y轴平行。z轴与图像平面垂直，为摄像机的光轴。z轴与图像平面的交点，O₀为图像坐标系的原点，由点O与x,y,z轴组成的直角坐标系称为摄像机坐标系。点O与O₀之间的距离为摄像机的焦距。X_w轴，Y_w轴，Z_w轴为在环境中选择的一个基准坐标系，即为世界坐标系。设空间中某点P在世界坐标系下的坐标是(X_w,Y_w,Z_w)，在摄像机坐标系下的坐标是(x_p,y_p,z_p)，像平面的理想坐标为(x_p′,y_p′)，受镜头畸变影响而偏移的实际像平面坐标为(x_p ^*,y_p ^*)，图像坐标为(u,v)。那么，从客观场景中的世界坐标到数字图像的成像变换需经历的步骤为：

步骤（1）：从客观场景的世界坐标(X_w,Y_w,Z_w)到摄像机坐标(x_p,y_p,z_p)的转换，该变换的齐次形式可表示为：

$\left[\begin{array}{c}{x}_p\\ y_p\\ z_p\\ 1\end{array}\right]=\left[\begin{array}{cc}R& T\\ 0& 1\end{array}\right]\left[\begin{array}{c}{X}_w\\ Y_w\\ Z_w\\ 1\end{array}\right],$

其中，R为3×3的旋转矩阵，该旋转矩阵由沿着x,y,z三个方向的轴的旋转所组成，若α，β，γ是旋转角度。则有

$R=\left[\begin{array}{ccc}\cos \left(\mathrm{\α}\right)& -\sin \left(\mathrm{\α}\right)& 0\\ \sin \left(\mathrm{\α}\right)& \cos \left(\mathrm{\α}\right)& 0\\ 1& 0& 1\end{array}\right]\left[\begin{array}{ccc}\cos \left(\mathrm{\β}\right)& 0& -\sin \left(\mathrm{\β}\right)\\ 0& 1& 0\\ \sin \left(\mathrm{\β}\right)& 0& \cos \left(\mathrm{\β}\right)\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \left(\mathrm{\γ}\right)& -\sin \left(\mathrm{\γ}\right)\\ 0& \sin \left(\mathrm{\γ}\right)& \cos \left(\mathrm{\γ}\right)\end{array}\right],$

T为平移向量。

步骤（2）：从摄像机坐标(x_p,y_p,z_p)到理想的像平面坐标(x_p′,y_p′)的转换，该变换可表示为如下：

${x_p}^{\′}=\frac{x_p}{z_p}f$

${y_p}^{\′}=\frac{y_p}{z_p}f$

则，用齐次矩阵的形式表示为：

$\left[\begin{array}{c}{x_p}^{\′}\\ {y_p}^{\′}\\ 1\end{array}\right]=\left[\begin{array}{cccc}\frac{f}{z_p}& 0& 0& 0\\ 0& \frac{f}{z_p}& 0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}{x}_p\\ y_p\\ z_p\\ 1\end{array}\right]$

步骤（3）：从理想的像平面坐标(x_p′,y_p′)到受镜头畸变影响而偏移的实际像平面坐标为(x_p ^*,y_p ^*)的转换，可用如下形式表示：

x_p ^*＝x_p′-φ_x ^*

y_p ^*-y_p′-φ_y ^*

其中，φ_x*，φ_y*为镜头畸变偏差。

步骤（4）：从图像坐标(u,v)与理想的像平面坐标(x_p',y_p')之间的变换。一幅图像可以被定义为一个二维函数f(u,v)，其中u和v是空间（平面）坐标，f在任何坐标点(u,v)处的振幅为图像在该点的亮度。假设图像的大小为M×N，在M行N列的图像中的每一个元素坐标(u,v)分别是该像素在数组中的列数与行数，则图像坐标系uv是以像素为单位的坐标系（如附图3示）。以图像内某一点O₀为原点，x'轴，y'轴分别与u轴，v轴平行，其中，(u,v)是图像坐标系中像素级单位的坐标，(x_p',y_p')是图像坐标系中毫米级（mm）单位的坐标。在x'y'坐标系中，原点O₀是摄像机光轴与图像平面的交点，设原点O₀在u,v坐标系中的坐标为(u₀，v₀)，dX与dY分别是每个像素在X轴与Y轴方向上的像素级的物理尺寸。而像平面坐标系x'y'是以物理单位（如mm）表示的坐标系。则当考虑摄像机畸变时，两坐标系之间的转换关系为：

$u=\frac{{x_p}^{\′}}{\mathrm{dX}}+{{\mathrm{\ηy}}_p}^{\′}+u_0$

$v=\frac{{y_p}^{\′}}{\mathrm{dY}}+v_0$

式中，η为像素点畸变的畸变因子，它是由于摄像机镜头光轴不垂直于CCD成像平面等因素造成的。

因此，当不考虑摄像机成像畸变时，客观场景的世界坐标、摄像机坐标、像平面坐标、图像坐标的关系可综合表示成下式：

$z_p\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{ccc}\frac{1}{\mathrm{dX}}& 0& u_0\\ 0& \frac{1}{\mathrm{dY}}& v_0\\ 0& 1& 0\end{array}\right]\left[\begin{array}{cccc}f& 0& 0& 0\\ 0& f& 0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cc}R& T\\ 0& 1\end{array}\right]\left[\begin{array}{c}{X}_w\\ Y_w\\ Z_w\\ 1\end{array}\right]$

$=\left[\begin{array}{cccc}\frac{f}{\mathrm{dX}}& 0& u_0& 0\\ 0& \frac{f}{\mathrm{dY}}& v_0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cc}R& T\\ 0& 1\end{array}\right]\left[\begin{array}{c}{X}_w\\ Y_w\\ Z_w\\ 1\end{array}\right]M_{\mathrm{internal}}{M}_{\mathrm{external}}{W}_{\mathrm{XYZ}}$

由上式可知，M_internal只与摄像机内部结构有关，其中x方向尺度因子y方向尺度因子摄像机光轴与图像平面的交点(u₀,v₀)，畸变因子η成为摄像机的内部参数；M_external由摄像机的外部参数决定，只与摄像机相对于世界坐标系的位置有关，主要包括旋转矩阵R的3个偏转角α，β，γ和平移向量T的3个参数，W_XYZ为客观场景中的点在世界坐标系下的齐次坐标。

因此，如果已知M_internal，M_external矩阵，对任何空间点，知道其世界坐标，就可以求出其图像像素坐标系中该图像点的坐标。但是，若只是已知某点的图像坐标，即便是M_internal，M_external矩阵已知，也不能唯一确定该点在客观场景中的世界坐标，因其该点对应客观场景空间中的一条射线。若考虑摄像机镜头的光学畸变因素，就要考虑下式：

$\left[\begin{array}{c}{x_p}^{\′}\\ {y_p}^{\′}\end{array}\right]=\left[\begin{array}{c}{x_p}^{*}\\ {y_p}^{*}\end{array}\right]+\left[\begin{array}{c}{{\mathrm{\φ}}_x}^{*}\\ {{\mathrm{\φ}}_y}^{*}\end{array}\right]$

根据以上所述的摄像机的透视投影模型，结合两高速摄像机的内外参数关系，运行matlab，可以计算出织物上各特征点的世界坐标[x,y,z]。数据如附图4～8所示，图4中第一行为其中的8个特征点在X方向上的坐标，第二行为其中的8个特征点在Y方向上的坐标，图5与图4类似，图6中第一、二、三行分别为X、Y、Z方向上其中的8个特征点的坐标。图7、图8分别为每个采样点在左相机、右相机图像中计算的可靠性，取值范围[-1，1]，数值越大，则计算结果越可靠。

第④步骤计算各特征点的基本运动量中，设任意特征点在一个运动周期内的第一帧的坐标为[x₁，y₁，z₁]，第i-1帧时的坐标为[x_i-1，y_i-1，z_i-1]，在第i帧时的坐标为[x_i，y_i，z_i]，在第i+1帧时的坐标为[x_i+1，y_i+1，z_i+1]则该特征点在第i-1、i帧之间的位移为 $s_i=\sqrt{{(x_i-x_{i-1})}^2+{(y_i-y_{i-1})}^2+{(z_i-z_{i-1})}^2},$ X方向上的位移为x_i-x_i-1，Y方向上的位移为y_i-y_i-1，Z方向上的位移为z_i-z_i-1。在第i、i+1帧间的位移为 $S_{i+1}=\sqrt{{(x_{i+1}-x_i)}^2+{(y_{i+1}-y_i)}^2+{(z_{i+1}-z_i)}^2},$ X方向上的位移为x_i+1-x_i，Y方向上的位移为y_i+1-y_i，Z方向上的位移为z_i+1-z_i。则到i帧时，Z方向上的累积位移为z_i-z₁，累积位移是指某特征点从图像序列开始运动到该图像序列时的总位移。在本实验中，拍摄频率为200帧/秒，因此两帧之间的时间间隔很小，为t＝0.005秒，则可把两帧之间的平均速度看做是后一帧时的即时速度，则第i帧时的速度为v_i＝s_i/t，第i+1帧时的速度为v_i+1＝s_i+1/t；III.根据相邻两帧的速度，可以求得从i帧到i+1帧间的加速度为从i+1帧到i+2帧间的加速度为则从加速度a_i+1变化到加速度a_i+2的加速度为 ${\mathrm{aa}}_{i+2}=\frac{a_{i+2}-a_{i+1}}{t}.$

第⑤步骤中提取了织物三个动态特征参数：每种织物能在Z方向上能达到的最大累积位移；“波”形特征曲线在最凹处时的曲率；“波”形曲线在运动过程中的保形率，即每一帧视图中“波”形曲线的相似度。附图10与附图11所示分别为采用附图1的拍摄装置针对涤纶纺丝软缎、真丝软缎、真丝绸、顺纡绉、平纹全棉、涤塔夫进行动态拍摄的部分运动周期（其中涤纶纺丝软缎的第81帧、真丝软缎的第37帧、真丝绸的606帧、顺纡绉的第31帧、平纹全棉的第30帧、涤塔夫的第68帧为“波”状曲线最明显时刻）、各织物织物的三种运动特征参数的计算结果。计算步骤如下：

i.Z方向上的最大累积位移

在每一幅图中，均存在一个特殊点，该点具有的特征为：X、Y、Z方向的合加速度（该加速度大于0）最小，速度增加到最大，则该点在Z方向的累积位移z最大。此时，在视图中，该点经过的“波”形曲线形状最明显。采用这样的方法，可以找到这样的一个特殊点。

ii.“波”状曲线拟合

涤纶仿丝软缎、真丝软缎、真丝绸、顺纡绉、平纹全棉、涤塔夫这六种织物所表现出的“波”状曲线有着明显的区别。涤纶仿丝软缎、真丝软缎的“波”状曲线比较圆润；真丝绸的“波”状曲线比较含蓄；顺纡绉、平纹全棉、涤塔夫的“波”状曲线比较生硬。为了定量地描述及区分这两种“波”状曲线，有必要对图像的“波”形做深入的研究，经定性分析，组成该条“波”形曲线的所有点运动到该图像序列时在Z方向上的累积位移最大。以某一帧为例，采用SPSS19.0对该视图上每个特征点的最大累积位移z_cd进行聚类，将那些有着较大累积位移的特征点提取出来。对应这些被选出来的特征点的X、Y坐标，“波”状曲线就在x-y平面上被拟合成二次曲线。就可以根据拟合出来的曲线算得每条“波”曲线在最低点时的曲率。

iii.“波”状曲线保形率

一般认为，织物的某一运动形态在传播过程中受织物这一介质的影响，如纱线间的紧密程度、纱线本身的物理机械性能。如当纱线间较为稀疏时，则运动形态的传播受纱线间的空气阻力的影响较大，在运动形态的传播过程中则表现为运动形态的保持差且很快消失。那么，在该试验中，研究“波”形曲线在各序列图像中的保形率（即各个序列图像中“波”形曲线的相似程度）就显得有意义了。曲线保形率也即各“波”状曲线在每个视图中的相似程度。根据拟合出的二次曲线表达式，在某一点的曲线曲率就可求得。曲率公式为假设某一幅视图中在某一点上的曲率为ρ_i(i＝1……n)，n为视图的序列号。因此，“波”状曲线的保形率为（其中那么，当保形率α越小时，在各视图中的曲线的相似程度越高。

步骤⑥为对织物的结构参数及物理性能进行测量。主要运用FAST织物风格仪、静态悬垂仪器、密度镜等，所测数据如附图17所示。其中W为织物的克重（单位：g/m²）、N_m为织物的纱支（单位：m/g）、P_J为织物的经向密度（单位：根/10cm）、P_W为织物的纬向密度（单位：根/10cm），E_J为织物的经向紧度（单位：%）、E_W为织物的纬向紧度（单位：%）、E_Z为织物的总紧度（单位：%），ST为织物的表观厚度（单位：mm）、B_J为经向弯曲刚度（单位：μN.m）、B_W为纬向弯曲刚度（单位：μN.m）、G为织物的剪切刚度、F_J为经向成形性（单位：mm²）、F_w为纬向成形性（单位：mm²）、δ为悬垂系数、ω为悬垂形态比系数、λ为悬垂形态系数比、为平均波峰半径、为平均波谷半径、为平均波峰间距、为平均波谷间距、Dσ_peaks为波峰间距方差、Dσ_valleys为波谷间距方差。

附图9为静态悬垂仪测出的投影图，其中1为波峰半径R_peak；2为波谷半径R_valley；3为波峰间距D_peaks；4为波谷间距D_valleys。

悬垂系数其中A_F为投影面积，A_d为圆盘面积，A_D为织物面积；

平均波峰半径n为波峰（或波谷）数，R_peak(i)为第i个波峰的半径；

平均波谷半径n为波谷（或波峰）数，R_valley(i)为第i个波谷的半径；

平均波峰间距n为波峰（或波谷）数，D_peak(i)为第i个波峰的间距；

平均波谷间距n为波谷（或波峰）数，D_valley(i)为第i个波谷的间距；

波峰间距方差n为波峰（或波谷）数，D_peak(i)为第i个波峰的半径,为平均波峰间距；

波谷半径方差n为波波谷（或波峰）数，D_valley(i)为第i个波谷的半径,为平均波谷间距；

第⑦步骤为基于织物结构参数及织物物理机械性能提取织物的形态风格的3个因子，利用SPSS实现，如附图18中所示，为因子分析后的解释的总方差表，前3个成分的初始特征值分别为11.020、4.179、3.776，均大于1，因此认为这3个因子得到的信息足以证明保留，而前三个成分特征值的累计贡献率为79.063%，因此，可以提取这3个因子。同时，如附图19所示，第一公因子更能代表悬垂系数δ、平均波峰间距平均波谷间距经向成形性F_J、经向弯曲刚度B_J、经向弯曲长度C_J、表面厚度ST、纬向弯曲刚度B_W、纱支N_m、总紧度E_Z、经向紧度E_J、纬向紧度E_W、纬向弯曲长度C_W，可以表现为织物形态的悬垂系数（主要为投影覆盖度）方面；第二公因子更适合代表纱支N_m、纬向紧度E_W、纬向弯曲长度C_W、经向密度P_J、克重W、纬向密度P_W、悬垂形态比系数ω，主要表现为悬垂的美感；第三公因子则更能代表纬向弯曲刚度B_W、总紧度E_Z、经向紧度E_J、纬向弯曲长度C_W、纬向密度P_W、纬向成形性F_W、剪切刚性G、波谷间距方差Dσ_valley、平均波峰半径波峰间距方差Dσ_peaks，主要表现为波的均匀、相似程度。进而，可以将这三个因子分别命名为：悬垂投影因子、悬垂美感因子、悬垂波状均匀度因子。则根据附图19中所示，各因子基于织物的结构参数、织物的物理机械性能的计算方程如下式(1)(2)(3)：

第⑧步骤为完成织物在视觉方面的评价，为了更有说服力，故提出的评价因子沿用基于织物结构参数、织物物理机械性能的评价方法提取出的评价因子，即步骤⑦中所述的评价因子。首先，建立每个织物的动态特征参数与织物物理性能间的关联关系，如附图20至22所示，然后根据（1）（2）（3）的3个因子式的得分方程式和关联关系，可得出基于三个动态特征参数（Z方向上的最大累积位移、“波”状曲线在明显时的曲率、曲线的保形率）的3个因子式的得分方程式，即将因子评价式中的每个织物物理性能参数都用运动特征参数的关系式来表示，如式（4）（5）（6）所示：

根据上述（4）（5）（6）三个方程式计算得到动态的体现织物的形态风格的3个因子，完成对织物的风格评价。

本发明分别基于织物静态物理性能和动态序列图像两方面对织物的风格进行评价，从织物形态的角度提出了三个在静态和动态都具体体现织物形态风格的可量化因子，分别是：悬垂投影因子、悬垂美感因子和波均匀度因子。在现在的单纯基于织物的物理性能进行外观形态评价时，需要对织物分指标多次测量，这种测量方法比较繁琐而且主观性太强，且基本为静态时测得；而在本发明的基于动态序列图像的织物风格评价方法中，只要对织物进行一次物理性能测量，然后从动态序列图像中分析计算得到各运动特征参数，建立织物的动态特征参数与织物静态物理性能间的关联关系，根据关联关系得到基于动态序列图像计算的动态特征参数的三个体现织物形态风格的主因子，即可完成对织物的风格评价。本发明解决了人眼视觉（即风格心理量）与静态物理量的联通关系，提出了三个无论在静态还是动态都能具体表征织物形态的特征参数（即三个形态风格评价主因子），实现了快速、简便、客观地完成对织物的风格评价。

Claims (6)

1.一种基于动态序列图像的织物形态风格评价方法，其特征在于：包括利用设置的多台CCD数码相机对织物面料的运动过程进行动态拍摄，得到织物运动的多帧动态序列视图的步骤；利用计算机视觉处理技术从同一帧的所有拍摄的视图上提取特征点的步骤；将对应视图上的特征点进行立体匹配并计算匹配点在世界坐标系中的三维坐标的步骤；从获得的特征点的三维坐标计算织物的基本运动物理量的步骤，所述基本运动物理量包括各特征点在X,Y,Z方向上的位移、速度、加速度；根据织物的所述基本运动物理量及运动特征提取织物动态特征参数的步骤，所述织物动态特征参数包括织物面料在Z方向上能达到的最大累积位移、运动过程中织物表现出来的“波”形曲线在最凹处时的曲率以及“波”形曲线在运动过程中的保形率；采用传统的测量方法对织物的物理性能进行测量的步骤，所述的物理性能包括织物结构参数与物理机械性能；然后，对这些测量的物理性能进行因子分析的步骤,提取计算与织物形态风格相关的3个主因子，所述3个主因子分别为悬垂投影因子、悬垂美感因子、波均匀度因子；计算所述动态特征参数与织物物理性能间的关联关系，结合对测量的物理性能进行因子分析的步骤中提取计算的形态风格评价因子，得到基于动态序列图像的3个形态风格评价主因子，完成织物形态在视觉方面的形态风格评价。

2.如权利要求1所述的基于动态序列图像的织物形态风格评价方法，其特征在于：所述传统的测量方法包括使用FAST织物风格仪、静态悬垂仪和密度镜对织物的物理性能进行测量。

3.如权利要求1所述的基于动态序列图像的织物形态风格评价方法，其特征在于：所述提取特征点的步骤中，采用Canny边缘检测算法和Harris角点检测算法结合进行特征点的提取，再采用SIFT特征点匹配算法和双目视觉视差原理获得各特征点在世界坐标系中的三维坐标。

4.如权利要求1所述的基于动态序列图像的织物形态风格评价方法，其特征在于：所述得到的基于动态序列图像的3个形态风格评价主因子，完成织物形态在视觉方面的评价，其中所述3个主因子的计算式分别如下：

其中，W为织物的克重(单位：g/m²)、N_m为织物的纱支(单位：m/g)、P_J为织物的经向密度(单位：根/10cm)、P_W为织物的纬向密度(单位：根/10cm)，E_J为织物的经向紧度(单位：％)、E_W为织物的纬向紧度(单位：％)、E_Z为织物的总紧度(单位：％)，ST为织物的表观厚度(单位：mm)、B_J为经向弯曲刚度(单位：μN.m)、B_W为纬向弯曲刚度(单位：μN.m)、C_J为经向弯曲长度，C_W为纬向弯曲长度，G为织物的剪切刚度、F_J为经向成形性(单位：mm²)、F_W为纬向成形性(单位：mm²)、δ为悬垂系数、ω为悬垂形态比系数、为平均波峰半径、为平均波谷半径、为平均波峰间距、为平均波谷间距、Dσ_peaks为波峰间距方差、Dσ_valleys为波谷间距方差，W1、Nm₁、ST₁、 G₁、δ₁、ω₁、 分别为各物理性能与最大累积位移Z的数量关系值,W₂、Nm₂、 ST₂、G₂、δ₂、ω₂、分别为各物理性能与波形曲线曲率的数量关系值,W₃、Nm₃、 ST₃、G₃、δ₃、ω₃、 分别为各物理性能与波形曲线保形率的数量关系值。

5.如权利要求1所述的基于动态序列图像的织物形态风格评价方法，其特征在于：将待评测的多种不同的织物面料进行裁剪，裁剪后将各织物面料的运动过程进行单个的多帧动态视图拍摄，提取出织物动态特征参数，依据参数进行对比和评价，所述待评测的多种不同的织物面料采用同一条件进行拍摄和织物动态特征提取。

6.如权利要求5所述的基于动态序列图像的织物形态风格评价方法，其特征在于：所述同一条件包括各织物面料平面尺寸的大小相同、触发面料运动的机制条件相同。