这周看了一篇Ryutarou Ohbuchi网格水印的论文，论文中提出在网格的频率域中加入水印。对于网格而言，没有如图像中的DCT等转换到频率域的变换，因此用什么量来模拟传统频率域中的系数，是很关键的。本文提出的水印算法在给定的3D多边形网格的变换域中嵌入水印，该变换域在本文中指网格频谱分析。同时，3D多边形网格的形状是由它的点连接关系和顶点坐标定义的。

文中通过更改网格频谱系数来实现水印的嵌入。先将原始网格变换到频率域中，然后修改频率域系数，嵌入水印信息，最后将模型从频率域反变换到空域即可得到嵌入水印的模型。

首先，介绍网格的频率域

计算网格频谱涉及到拉普拉斯矩阵特征值分解。从特征值分解中会产生一个特征值序列和对应的特征向量序列。小的特征值对应空间频率中的低频，大的特征值对应空间频率中的高频。小的特征值对应的特征向量和低频系数代表了全局形状特征，而大的特征值对应的特征向量和高频系数代表局部或细节形状特征。将模型中的空域中的点投影到归一化的特征向量将生成该顶点的网格频率系数。文中使用联合拉普拉斯矩阵，即基尔霍夫矩阵K来替代原始的拉普拉斯矩阵。

基尔霍夫矩阵K定义如下：（实际上这个矩阵就是只用到顶点连接关系的拉普拉斯矩阵）

K = D – A

D：对角矩阵，它的元素D_ii = d_i，表示节点i的度数，

A：多边形网格的邻接矩阵，其元素a_ij定义如下：

（1）

一个有n个顶点的多边形网格M将会产生一个n x n大小的K矩阵，其特征值分解将产生n个特征值λ_i 和 n 个n维特征向量 w_i，（1 ≤ i ≤ n ）。将原始网格中的第i个顶点的每个分量 Vi = （x_i， y_i， z_i）（1 ≤ i ≤ n ）单独地投影到第i个归一化的特征向量e_i上，即e_i = w_i / ||w_i||。（也就是说e_i是w_i的单位向量）。

通过以上步骤产生了n个网格频率系数 ri = ( rs,i ,rt,i ,ru,i ) ，（1 ≤ i ≤ n ）。下标s，t和u表示的是在网格频率域中的正交坐标轴，对应空间坐标轴x，y和z。为了从频率域中逆变换到空间域，需要将e_i与r_i相乘并求和，具体转换公式如下：

（2）

水印生成

上文中已经提到过，文中通过更改网格频谱系数来实现水印的嵌入。每个频率系数有三个频率轴s，t和u，通过对频率系数的三个分量进行调制以完成对水印的嵌入。

本文嵌入的数据是一个m维的比特向量a = (a1,a2,…,am) , 。每一个比特aj乘以码片速率c得到水印信号向量b = (b1,b2,…,bmc) , ，且b的总长度为mc

bi = aj , j.c ≤ i < (j+1).c ( 3 )

（注：我对这个转换的理解是，将原本要嵌入的序列的每一位都重复嵌入c（码片速率）次，比如，m = 5，c = 3，原始序列为a =（11001），经过上式转换将变为 b =（111，111，000，000，111），文中也解释说将同一个比特重复嵌入c次可以增加水印抵抗随机调价噪声攻击）

将b转换为b'， ，转换规则为

（4）

水印嵌入

接下来就可以考虑嵌入水印的操作了。首先以频率域中的s轴为例进行水印嵌入。假定r_si是嵌入水印前的网格的第i个频谱系数的s轴分量，pi∈{1,-1}是由密钥Kw生成的一个伪随机序列，α (α> 0)是调制幅度。则第i个嵌入水印后的频谱系数计算公式如下：

（5）

在水印提取过程中，需要使用水印嵌入过程中用到的密钥，也就是作为伪随机数的种子值要相同。因此，密钥的分发可以通过如公钥加密等分方式来完成。

在频谱的t轴和u轴都执行如上的操作后将生成一系列嵌入水印的频谱系数，使用上文中的公式（2）将频谱系数做逆变换即可得到空域中的顶点坐标，如此便可得到一个嵌入水印的多边形网格。

水印提取

作为一个私有水印，本文中提出的水印算法是非盲水印算法，也就是说，在水印提取阶段需要原始的未嵌入水印的网格，记为cover-mesh、和嵌入水印的（或者被降级的）网格，记为stego-mesh。水印提取时，首先对cover-mesh M和stego-mesh Mˆ，进行对齐。（嵌入水印后的三维模型在使用的时候经常要通过一定的相似变换处理，如网格模型的旋转、平移、缩放等操作。这些操作带来了模型顶点信息的改变，而三维模型数字水印算法中水印信息就是嵌入在顶点的属性信息中，在此情况下是不能提取出嵌入的水印信息的。为此，在提取水印信息之前需要把待检测模型变换到原始网格模型的坐标下。这一步骤称为网格对齐操作-----来自硕士论文《三维模型数字水印技术研究》）

为了对齐两个网格，文中使用网格前5个（最低频）频谱系数来粗略地重构出网格的形状。然后使用重建的形状得到一个3 x 3 的协方差矩阵，由该矩阵计算出一系列的特征向量。通过比较由两个形状计算出的特征向量将两个网格对齐。

对齐后的网格各自进行特征值分解，以得到各自的频谱系数。其中，M分解得到rsi，M^分解得到rsi^,计算。

计算三个频率轴上的系数总和qj

（6）

上式中使用的p_i与嵌入时保持一致。同时，假定M^中的顶点坐标受到的扰动可以忽略，那么

（7）

其中，q_j的取值为{αc , -αc }.由于α和C始终为正，简单测试q_j的符号即可恢复原始的信息比特序列a_j

a_j = sign(q_j)

上式中的sign为取符号函数。通过本文中的公式（4）可以转换得到原始的信息比特序列b_i。

优化：网格分割

特征值分解对于只有几百个顶点的网格来说效果较好，但是当网格规模变大时，如顶点个数为10^4 ~ 10^7时，先前的方法就会出现两个问题：一是计算特征值分解的时间远大于水印嵌入时间，二是分解数值的稳定性变得越来越有问题.

因此,当要对顶点数目为10^4 ~ 10^7的网格进行操作时,需要将原始网格分割成小的子网格(如500个顶点),从而水印的嵌入和提取都在每个子网格中独立地进行。

网格分割有一个额外的好处。    将水印信息重复地嵌入到子网格中，可以抵抗网格切除等攻击，只要被切除的网格中保留了至少一个之前划分的子网格，水印信息就可以从中提取出。

本文实现了一个简单的网格分割。首先，手动地选取"特征点"，使得特征点近似在预期的子网格的中心位置，并使特征点近似均匀分布在整个网格上。每个子网格沿着特征点根据网格顶点之间的连接关系逐步进行区域的扩大。当原始网格完全被分割时，子网格停止增长。

网格频谱分析在各个子网格中进行，忽略子网格边界上顶点的连接。在网格的嵌入和提取阶段要求子网格划分情况相同，因此在原始网格中用到的特征点需要保存。

其他补充

由于在水印提取阶段用到了最低频的5个频率系数进行网格的对齐，因此在嵌入水印的时候，并不是将水印嵌入到所有的频率系数中的，最多只能嵌入到 n – 5 个高频系数中。