CN101833742A - 基于球面参数化的三维网格模型数字水印方法 - Google Patents

Abstract

本发明涉及一种基于球面参数化的三维网格模型数字水印方法。当前数字水印技术大都是针对静止图像、视频流和音频流这些媒体数据类型的，而对三维几何模型数据的水印技术的研究工作相对较少。本发明的三维网格模型水印方法利用网格全局球面参数化的方法得到网格的球面调和分析系数，将水印嵌入调和分析系数中，然后进行相应的逆变换得到嵌入水印的网格。本方法还解决了由于采样和嵌入水印而造成的网格变形问题，而且在提取水印过程中将网格水印信息进行放大和对水印的特殊编码策略以利于水印的提取。本发明方法提取水印不需要原始网格，不需要预处理，可以抵抗很多攻击，攻击包括加噪、光滑、增强、旋转、平移、缩放、重采样、裁剪和一些组合攻击。

Description

基于球面参数化的三维网格模型数字水印方法

技术领域

本发明涉及一种基于球面参数化的三维网格模型数字水印方法。

背景技术

当前数字水印技术大都是针对静止图像、视频流和音频流这些媒体数据类型的，而对三维几何模型数据的水印技术的研究工作相对较少。随着三维扫描技术的发展，三维几何模型已经成为继音图视频后的一种新的媒体类型，被广泛应用在娱乐业和制造工业以及其他各种领域中。

由于三维模型数据自身的特点，使得传统的图像水印方法不能简单照搬地应用于三维几何模型，三维模型水印存在一些难点，比如由于三维几何模型数据具有不规则性，所以在水印嵌入过程中，缺乏进行频率分解的某种自然地参数化方法。三维模型数据由点、线、面等要素构成，这些要素可以组合成各种不同的数据表达方式，另外三维模型数据的各要素集合没有一个固定的排序标准，而对静止图像来说却可以按照像素点的平面位置排序，对音频流和视频流数据可以按照时间轴来排序，对于这种不规则的数据类型，不能简单地应用已有的各种变换域水印方法，需要寻找适当的能够反映三维模型数据特征的参数用于各种变换域水印方法。另外，水印检测过程中，几何模型简化操作和其它的攻击方法可能会改变几何模型的连接关系，或拓扑连接。

与静止图像、音频、视频相比，三维模型的表示具有如下特点：

1.没有固有的数据顺序。音频、视频的数据是按时间顺序排列的，静止图像、视频的帧则以扫描线顺序排列，而三维模型却没有一个固定的数据顺序。

2.没有明确的采样率的概念，三维表面模型中的数据，不具有像图像、音频、视频那样的方便的数学工具(如余弦变换、Fourier变换、小波变换等)可以使用。

3.模型数据中不但包含几何信息还有拓扑信息，使水印提取时候的同步问题更加复杂。

4.表示方法不唯一。同一个三维模型，可以用多种不同的模型表示，在不同模型间进行转换的过程中，容易引起属性损失。没有图像、视频中自然存在适合于嵌入水印的噪声存在。

5.对模型进行几何和拓扑操作的工具很多，使对模型的修改更加容易。一些图像水印、音频水印、视频水印中的问题，如有损压缩、同步问题等。三维水印同样存在。

表示三维模型的主要方法有：

线框模型：它是利用对象形体的棱边和顶点来表示几何形状的一种模型，只能反映出三维物体的一部分形状信息，难以表示物体的剖面图、消隐线及两个面的交线或轮廓线等，不能完整的定义三维物体。线框模型在三维方面的处理上有很多困难，如消隐、着色、特征处理等。

表面模型：它把线框模型中棱线所包围的部分定义成形体的表面，增加面的有关信息及连接指针，然后利用形体表面的集合来描述形体的形状。表面模型可分为多边形网格表示法和参数曲面表示法。多边形可以使用足够的细节，可以表示任何表面和描述任意复杂的物体。事实上，没有多少物体不能使用多边形网格表示。而参数曲面表示法用少量的控制顶点及形状参数来描述物体表面的形状。常见的有Bezier、B样条等参数曲线、参数曲面造型技术及元球(Metaball)隐式曲面造型技术。参数曲面与多边形网格相比，更简洁，局部修改容易，精度控制方便。但对于复杂的局部表面区域，参数曲面的拼接往往因涉及到复杂的约束求解而使人望而却步。目前，已有一些方法成功地将多边形网格局部参数化，可以方便地处理多边形网格表面。

实体模型：在这种表示方法中，把许多具有一定形状和体积的基本体素(如立方体、圆柱、圆锥等)通过并、交、补等布尔运算和基本变形操作建立的三维立体模型，实体建模的特点在于建立一个物体的完整的形状模型，它有明确的物体包容空间，且各表面间有严格的拓扑关系，形成一个整体。实体模型技术只能表示简单的三维物体，而不能表示复杂物体的表面。实体模型理论的发展可以追溯到1970年，当时是利用CSG(Constructive Solid Geometry)方法。CSG建模方法其实是将最基本的实体(立方体、圆柱体、圆锥体等)进行布尔运算，这就需要事先按一定的顺序建立好大小、位置合适的基本实体，并且不能改变。但实际情况是谁也不能保证设计结果不被修改，所以CSG方法不能被设计人员接受。

三维体数据表示法：现实世界中的许多物体、自然现象及一些计算模型，其表面和内部都包含着对象信息，不能用曲线或曲面等几何造型方法表示，而只能用三维体数据表示。

同图像水印方法分类相似，三维网格模型水印方法也可以大致分为两类，一类是空间域水印方法，一类是变换域水印方法。早先的三维模型水印方法研究也是从空间域水印方法开始的，现在逐步向变换域水印方向发展。空间域网格水印方法对噪声和裁剪等攻击普遍没有抵抗力。

发明内容

本发明的目的在于提供一种基于球面参数化的三维网格模型数字水印方法。

本发明的三维网格模型水印方法利用网格全局球面参数化的方法得到网格的球面调和分析系数，将水印嵌入调和分析系数中，然后进行相应的逆变换得到嵌入水印的网格。本方法还解决了由于采样和嵌入水印而造成的网格变形问题，而且在提取水印过程中将网格水印信息进行放大和对水印的特殊编码策略以利于水印的提取。

本发明方法的具体步骤是：

步骤(1)对原始网格X-坐标进行重新排序，使其成为由小到大上升排列方式得到一个新X-坐标的网格；

步骤(2)采用球面参数化方法将原始网格M上的几何位置(即顶点坐标)信号F_t转换成球面参数网格上的信号F_t ^s(θ，φ)，t＝1，2，3，t表示三个坐标方向；θ(π≤θ≤0)表示球面上的经度角、φ(2π≥φ≥0)表示球面上的纬度角；

步骤(3)将F_t ^s(θ，φ)在经纬θ和φ方向做均匀采样得到规则采样信号F′_t，采样率为64×64。其中采样点θ_j＝π(2j+1)/4B，φ_k＝2πk/2B，B＝32，j、k为采样点序列号；

步骤(4)对F′_t做球面调和分析，得到F′_t的调和系数I_t，设单位球面上任意点的位置p用经度角θ和纬度角φ确定，即p＝(cosφsinθ，sinφsinθ，cosθ)，任意一球面信号表示为二维连续函数f(θ，φ)。球面调和分析的基函数(也被称为球面调和函数)Y_l ^m(θ，φ)表示为：

$Y_l^m(\mathrm{\θ},\mathrm{\φ})=K_{l,m}{P}_l^m{\left(\cos \mathrm{\θ}\right)e}^{\mathrm{im\φ}}---\left(1\right)$

m＝-l，-l+1，...，l-1，l；l＝0，1，2，...，

其中

是规一化常数，P_l ^m(x)是度为l、次为m的联合Legendre多项式，e为指数、i为复数，球面调和函数来自球面上的Laplacian方程，它们构成了球面连续函数空间中的一组完备的单位正交基，即：

$<Y_l^m,Y_k^n>={\mathrm{\&delta;}}_{l,k}{\mathrm{\&delta;}}_{m,n}$ 其中， ${\mathrm{\&delta;}}_{l,k}=\left\{\begin{array}{cc}1& l=k\\ 0& l\&NotEqual;k\end{array}\right.,$ ${\mathrm{\&delta;}}_{m,n}=\left\{\begin{array}{cc}1& m=n\\ 0& m\&NotEqual;n\end{array}\right.$

δ为脉冲函数

因此，任意定义在球面上的连续函数f(θ，φ)都可以用球面调和函数展开：

$f(\mathrm{\&theta;},\mathrm{\&phi;})=\underset{l=0}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\underset{m=-l}{\overset{l}{\mathrm{\&Sigma;}}}\hat{f}(l,m)Y_l^m(\mathrm{\&theta;},\mathrm{\&phi;})---\left(2\right)$

其中

为f(θ，φ)和Y_l ^m(θ，φ)的内积，称为(l，m)-Fourier系数：

$\hat{f}(l,m)=<f,Y_l^m>=K_{l,m}{\&Integral;}_0^{\mathrm{\&pi;}}\left[{\&Integral;}_0^{2\mathrm{\&pi;}}{e}^{-\mathrm{im\&phi;}}f(\mathrm{\&theta;},\mathrm{\&phi;})\mathrm{d\&phi;}\right]\&times;P_l^m\left(\cos \mathrm{\&theta;}\right)\sin \mathrm{\&theta;d\&theta;}---\left(3\right)$

公式(3)就是球面调和分析公式，所得的

即为球面调和系数；设f(θ，φ)的带宽为B，公式3中的积分运算简化为2B×2B的采样数据的加权和：

$\hat{f}(l,m)=\frac{\sqrt{2\mathrm{\&pi;}}}{2B}\underset{j=0}{\overset{2B-1}{\mathrm{\&Sigma;}}}\underset{k=0}{\overset{2B-1}{\mathrm{\&Sigma;}}}{a}_j^{\left(B\right)}f({\mathrm{\&theta;}}_j,{\mathrm{\&phi;}}_k)e^{-\mathrm{im}{\mathrm{\&phi;}}_k}{P}_l^m\left(\cos {\mathrm{\&theta;}}_j\right)---\left(4\right)$

其中采样点θ_j＝π(2j+1)/4B，φ_k＝2πk/2B，a_j ^(2B)为权值；

步骤(5)将水印进行编码，就是将字符串w作为水印信息，字符串包括四种字符，每种字符是由5个0和1个1组成的一个序列，字符串w的比特流表示为

w＝(w₁，w₂，...，w₂₄)，w_i∈{0.1}，i＝1，2，...，24。

步骤(6)取球面调和系数I₁中某些固定位置进行嵌入水印，按照以下水印嵌入公式加入水印信息，得到加入水印序列W的调和系数I′₁，

I′₁(k，l)＝I₁(k，l)+θ        (5)

其中k和l是在I₁中水印嵌入位置的行数和列数，

k＝22，23，24，l＝34，35，…，42，w_i＝0时θ＝0，w_i＝1时θ＝1；

步骤(7)对I′₁，I₂和I₃进行逆调和分析，利用变换公式(2)得到嵌入水印的采样信号F_t ^*；

步骤(8)对F_t ^*做反采样和反球面参数化得到加入了水印的模型M′；反采样和反球面参数化方法利用步骤(2)和(3)相应的反变换方法；

步骤(9)调整嵌入水印后的网格

通过将嵌入水印的网格进行调整来减少变形，设原始网格为M，嵌入水印后的网格为M^#，将原始网格进行球面参数化、采样、反采样和反球面参数化而不嵌入水印得到的网格为M^*，调整后的嵌入水印的网格为M^。修正公式如下：

M^＝M+α×[M^#+(M^*-M)-M]，α＜1(6)   α为缩小系数；

步骤(10)提取嵌入的水印采用如下方法：

步骤a.对待检测网格M^*的X-坐标进行重新排序，使其成为由小到大上升排列方式得到一个新X-坐标的网格；

步骤b.将待检测网格M^*的X-坐标按照下式进行处理，得到网格M′：

其中，f(x)为待检测网格M^*的原始X-坐标，g(x)为处理后的待检测网格M′的原始X-坐标，T为设定的阈值，f(m)为f(x)以及前后各三个点的和

f(m)＝f(x-3)+f(x-2)+f(x-1)+f(x)+f(x+1)+f(x+2)+f(x+3)

将g(x)进行放大处理得到放大坐标f′(x)，放大函数定义如下：

f′(x)＝g(x)+(f(x)-g(x))×β      (8)

其中β为放大系数。

步骤c.利用原始网格的球面参数化的信息将网格M′的几何位置(即顶点坐标)信号F_t进行均匀采样得到采样信号F′_t，采样率为64×64，t＝1，2，3，t表示三个坐标方向；

步骤d.对F′_t做球面调和分析，得到F′_t的调和系数I_t，球面调和分析方法与嵌入水印的球面调和分析方法相同。

步骤e.将I₁中频那些嵌入水印的数值取出组成水印信息，就把水印恢复出来，即将球面调和系数中嵌入水印位置上的值取出组成序列W，将W分成四组，每组长度为6，将每组中5个最小值改为0，1个最小值改为1，则得到的编码根据水印的编码将水印信息恢复出来。

本发明方法提取水印不需要原始网格，不需要预处理。可以抵抗很多攻击，攻击包括加噪、光滑、增强、旋转、平移、缩放、重采样、裁剪和一些组合攻击。

具体实施方式：

基于球面参数化的三维网格模型数字水印方法的具体步骤是：

步骤(1)对原始网格X-坐标进行重新排序，使其成为由小到大上升排列方式得到一个新X-坐标的网格；

步骤(2)采用球面参数化方法将原始网格M上的几何位置(即顶点坐标)信号F_t转换成球面参数网格上的信号F_t ^s(θ，φ)，t＝1，2，3，t表示三个坐标方向；θ(π≤θ≤0)表示球面上的经度角、φ(2π≥φ≥0)表示球面上的纬度角。

上述的球面参数化方法采用现有的成熟方法，原理如下：

球面参数化的方法针对亏格非零的网格模型用单位球面来代替复杂的网格表面作为几何信号的定义域，这样适用于球面的一些正交分析工具(如球面调和分析)就可以被用来处理几何信号。球面参数化处理几何信号的流程如下：

通过为任意网格构造一个拓扑同构的球面网格，定义在原网格表面的信号就被转换为定义在球面网格上的球面信号。

球面参数化方法分为两个步骤：

①生成带有局部参数化信息的累进网格表示

循环地执行边折叠简化操作，直到当前简化网格变成一个凸多面体(这个凸多面体总是存在的，最坏的情况就是四面体)，这个凸多面体被称为基网格。对每次边折叠简化操作，被删除的两个顶点被分别局部参数化到该边折叠操作生成的简化网格表面上。这些局部参数化信息被记录在累进网格表示的相应顶点分裂操作中。

②从基网格的中心投影可以得到相应的球面网格

从初始球面网格开始，以逆序累进地执行记录在累进网格表示中的顶点分裂操作。对每次顶点分裂操作，使用局部参数化信息把两个分裂出来的顶点放置在单位球面上。当所有的顶点分裂操作执行完毕，原始网格的球面参数化表示就生成了。

步骤(3)将F_t ^s(θ，φ)在经纬θ和φ方向做均匀采样得到规则采样信号F′_t，采样率为64×64。其中采样点θ_j＝π(2j+1)/4B，φ_k＝2πk/2B，B＝32，j、k为采样点序列号。

步骤(4)对F′_t做球面调和分析，得到F′_t的调和系数I_t，设单位球面上任意点的位置p用经度角θ和纬度角φ确定，即p＝(cosφsinθ，sinφsinθ，cosθ)，任意一球面信号表示为二维连续函数f(θ，φ)。球面调和分析的基函数(也被称为球面调和函数)Y_l ^m(θ，φ)表示为：

$Y_l^m(\mathrm{\&theta;},\mathrm{\&phi;})=K_{l,m}{P}_l^m{\left(\cos \mathrm{\&theta;}\right)e}^{\mathrm{im\&phi;}}---\left(1\right)$

m＝-l，-l+1，...，l-1，l；l＝0，1，2，...，

其中

是规一化常数，P_l ^m(x)是度为l、次为m的联合Legendre多项式，e为指数、i为复数，球面调和函数来自球面上的Laplacian方程，它们构成了球面连续函数空间中的一组完备的单位正交基，即：

$<Y_l^m,Y_k^n>={\mathrm{\&delta;}}_{l,k}{\mathrm{\&delta;}}_{m,n}$ 其中， ${\mathrm{\&delta;}}_{l,k}=\left\{\begin{array}{cc}1& l=k\\ 0& l\&NotEqual;k\end{array}\right.,$ ${\mathrm{\&delta;}}_{m,n}=\left\{\begin{array}{cc}1& m=n\\ 0& m\&NotEqual;n\end{array}\right.$

δ为脉冲函数

因此，任意定义在球面上的连续函数f(θ，φ)都可以用球面调和函数展开：

$f(\mathrm{\&theta;},\mathrm{\&phi;})=\underset{l=0}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\underset{m=-l}{\overset{l}{\mathrm{\&Sigma;}}}\hat{f}(l,m)Y_l^m(\mathrm{\&theta;},\mathrm{\&phi;})---\left(2\right)$

其中

为f(θ，φ)和Y_l ^m(θ，φ)的内积，称为(l，m)-Fourier系数：

$\hat{f}(l,m)=<f,Y_l^m>=K_{l,m}{\&Integral;}_0^{\mathrm{\&pi;}}\left[{\&Integral;}_0^{2\mathrm{\&pi;}}{e}^{-\mathrm{im\&phi;}}f(\mathrm{\&theta;},\mathrm{\&phi;})\mathrm{d\&phi;}\right]\&times;P_l^m\left(\cos \mathrm{\&theta;}\right)\sin \mathrm{\&theta;d\&theta;}---\left(3\right)$

公式(3)就是球面调和分析公式，所得的

即为球面调和系数。

如果存在着某个正数B＞0使得对所有的l＞B都有则称f(θ，φ)是带宽有限的函数，或者说f(θ，φ)的带宽为B。由Nyquist采样定律，对带宽为B的二维信号，只需要2B×2B的采样数据就可以完全恢复原信号。因此，设f(θ，φ)的带宽为B，公式3中的积分运算简化为2B×2B的采样数据的加权和：

$\hat{f}(l,m)=\frac{\sqrt{2\mathrm{\&pi;}}}{2B}\underset{j=0}{\overset{2B-1}{\mathrm{\&Sigma;}}}\underset{k=0}{\overset{2B-1}{\mathrm{\&Sigma;}}}{a}_j^{\left(B\right)}f({\mathrm{\&theta;}}_j,{\mathrm{\&phi;}}_k)e^{-\mathrm{im}{\mathrm{\&phi;}}_k}{P}_l^m\left({\cos \mathrm{\&theta;}}_j\right)---\left(4\right)$

其中采样点θ_j＝π(2j+1)/4B，φ_k＝2πk/2B，a_j ^(2B)为权值，其作用类似于积分运算中的sinθ，公式4就是球面函数f(θ，φ)的离散Fourier变换。

球面调和分析保证本文水印嵌入在真正Fourier意义上的频率，对噪声、裁剪、滤波和增强的攻击有很强的强壮性；对噪声、裁剪、理想滤波、增强、相似变换、几何变换、网格重采样和一些组合攻击后的网格进行水印提取时不需要网格对准和重采样。

步骤(5)将水印进行编码，就是将字符串w作为水印信息，字符串包括四种字符，每种字符是由5个0和1个1组成的一个序列，字符串w的比特流表示为

w＝(w₁，w₂，...，w₂₄)，w_i∈{0.1}，i＝1，2，...，24。

为了更好地提取球面调和系数中的水印信息，提取水印时不需要提取参数。以嵌入水印是字符串为例，将每个字符转成长度为n的由0和1组成的编码，每个编码中有k个0。在提取水印时，将球面调和系数中嵌入水印的那些位置上的值取出组成序列W，将W分成L/n组，L为序列W总长度，每组长度为n，将每组中k个最小值改为0，其他改为1，则得到的编码根据水印的编码将水印信息恢复出来。

例如，在水印编码中，字符“a”的编码为[1 0 0 0 0 0]；字符“b”的编码为[0 0 0 0 1 0]；字符“c”的编码为[0 0 0 1 0 0]。若在球面调和系数中相应的位置上提取的数值为[123 45 67 89 200 189 230 113 77 92 11967 23 98 134 250 54 23]，则将其分为三组，每组长度为6，进一步将这三组分别转化为[0 0 0 0 1 0]、[1 0 0 0 0 0]、[0 0 0 1 0 0]，则根据编码得到恢复后的水印信息为“bac”。

步骤(6)取球面调和系数I₁中某些固定位置进行嵌入水印，按照以下水印嵌入公式加入水印信息，得到加入水印序列W的调和系数I′₁，

I′₁(k，l)＝I₁(k，l)+θ        (5)

其中k和l是在I₁中水印嵌入位置的行数和列数，

k＝22，23，24，l＝34，35，…，42，w_i＝0时θ＝0，w_i＝1时θ＝1；

步骤(7)对I′₁，I₂和I₃进行逆调和分析，利用变换公式(2)得到嵌入水印的采样信号F_t ^*；

步骤(8)对F_t ^*做反采样和反球面参数化得到加入了水印的模型M′；反采样和反球面参数化方法利用步骤(2)和(3)相应的反变换方法；

步骤(9)调整嵌入水印后的网格

嵌入水印后的网格变形由两方面原因引起，一方面，由于采样和反采样是用插值来实现的，故存在误差。另一方面，嵌入水印也会给网格带来比较大的变形。我们通过将嵌入水印的网格进行调整来减少变形，设原始网格为M，嵌入水印后的网格为M^#，将原始网格进行球面参数化、采样、反采样和反球面参数化而不嵌入水印得到的网格为M^*，调整后的嵌入水印的网格为M^。修正公式如下：

M^＝M+α×[M^#+(M^*-M)-M]，α＜1(6)        α为缩小系数；

步骤(10)提取嵌入的水印采用如下方法：

步骤a.对待检测网格M^*的X-坐标进行重新排序，使其成为由小到大上升排列方式得到一个新X-坐标的网格；

步骤b.将待检测网格M^*的X-坐标按照下式进行处理，得到网格M′：

其中，f(x)为待检测网格M^*的原始X-坐标，g(x)为处理后的待检测网格M′的原始X-坐标，T为设定的阈值，f(m)为f(x)以及前后各三个点的和

f(m)＝f(x-3)+f(x-2)+f(x-1)+f(x)+f(x+1)+f(x+2)+f(x+3)

将g(x)进行放大处理得到放大坐标f′(x)，放大函数定义如下：

f′(x)＝g(x)+(f(x)-g(x))×β        (8)

其中β为放大系数。

步骤c.利用原始网格的球面参数化的信息将网格M′的几何位置(即顶点坐标)信号F_t进行均匀采样得到采样信号F′_t，采样率为64×64，t＝1，2，3，t表示三个坐标方向；

步骤d.对F′_t做球面调和分析，得到F′_t的调和系数I_t，球面调和分析方法与嵌入水印的球面调和分析方法相同。

步骤e.将I₁中频那些嵌入水印的数值取出组成水印信息，就把水印恢复出来，即将球面调和系数中嵌入水印位置上的值取出组成序列W，将W分成四组，每组长度为6，将每组中5个最小值改为0，1个最小值改为1，则得到的编码根据水印的编码将水印信息恢复出来。

Claims (1)

1.基于球面参数化的三维网格模型数字水印方法，其特征在于该方法的具体步骤是：

方法的具体步骤是：

步骤(1)对原始网格X-坐标进行重新排序，使其成为由小到大上升排列方式得到一个新X-坐标的网格；

步骤(2)采用球面参数化方法将原始网格M上的几何位置信号F_t转换成球面参数网格上的信号F_t ^s(θ，φ)，t＝1，2，3，t表示三个坐标方向；θ(π≤θ≤0)表示球面上的经度角、φ(2π≥φ≥0)表示球面上的纬度角；

步骤(3)将F_t ^s(θ，φ)在经纬θ和φ方向做均匀采样得到规则采样信号F_t′，采样率为64×64；其中采样点θ_j＝π(2j+1)/4B，φ_k＝2πk/2B，B＝32，j、k为采样点序列号；

步骤(4)对F_t′做球面调和分析，得到F_t′的调和系数I_t，设单位球面上任意点的位置p用经度角θ和纬度角φ确定，即p＝(cosφsinθ，sinφsinθ，cosθ)，任意一球面信号表示为二维连续函数f(θ，φ)；球面调和分析的基函数Y_l ^m(θ，φ)表示为：

Y_l ^m(θ，φ)＝K_l，mP_l ^m(cosθ)e^imφ    (1)

m＝-l，-l+1，...，l-1，l；l＝0，1，2，...，

其中

是规一化常数，P_l ^m(x)是度为l、次为m的联合Legendre多项式，e为指数、i为复数，球面调和函数来自球面上的Laplacian方程，它们构成了球面连续函数空间中的一组完备的单位正交基，即：

$<Y_l^m,Y_k^n>={\mathrm{\&delta;}}_{l,k}{\mathrm{\&delta;}}_{m,n}$ 其中， ${\mathrm{\&delta;}}_{l,k}=\left\{\begin{array}{cc}1& l=k\\ 0& l\&NotEqual;k\end{array}\right.,$ ${\mathrm{\&delta;}}_{m,n}=\left\{\begin{array}{cc}1& m=n\\ 0& m\&NotEqual;n\end{array}\right.$

δ为脉冲函数

任意定义在球面上的连续函数f(θ，φ)用球面调和函数展开：

$f(\mathrm{\&theta;},\mathrm{\&phi;})=\underset{l=0}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\underset{m=-l}{\overset{l}{\mathrm{\&Sigma;}}}\hat{f}(l,m)Y_l^m(\mathrm{\&theta;},\mathrm{\&phi;})---\left(2\right)$

其中

为f(θ，φ)和Y_l ^m(θ，φ)的内积，称为(l，m)-Fourier系数：

$\hat{f}(l,m)=<f,Y_l^m>=K_{l,m}{\&Integral;}_0^{\mathrm{\&pi;}}\left[{\&Integral;}_0^{2\mathrm{\&pi;}}{e}^{\mathrm{im\&phi;}}f(\mathrm{\&theta;},\mathrm{\&phi;})\mathrm{d\&phi;}\right]\&times;P_l^m\left(\cos \mathrm{\&theta;}\right)\sin \mathrm{\&theta;d\&theta;}---\left(3\right)$

公式(3)就是球面调和分析公式，所得的

即为球面调和系数；设f(θ，φ)的带宽为B，公式(3)中的积分运算简化为2B×2B的采样数据的加权和：

$\hat{f}(l,m)=\frac{\sqrt{2\mathrm{\&pi;}}}{2B}\underset{j=0}{\overset{2B-1}{\mathrm{\&Sigma;}}}\underset{k=0}{\overset{2B-1}{\mathrm{\&Sigma;}}}{a}_j^{\left(B\right)}f({\mathrm{\&theta;}}_j,{\mathrm{\&phi;}}_k)e^{-\mathrm{im}{\mathrm{\&phi;}}_k}{P}_l^m\left(\cos {\mathrm{\&theta;}}_j\right)---\left(4\right)$

其中采样点θ_j＝π(2j+1)/4B，φ_k＝2πk/2B，a_j ^(2B)为权值；

步骤(5)将水印进行编码，就是将字符串w作为水印信息，字符串包括四种字符，每种字符是由5个0和1个1组成的一个序列，字符串w的比特流表示为

w＝(w₁，w₂，...，w₂₄)，w_i∈{0.1}，i＝1，2，...，24；

步骤(6)取球面调和系数I₁中某些固定位置进行嵌入水印，按照以下水印嵌入公式加入水印信息，得到加入水印序列W的调和系数I′₁，

I′₁(k，l)＝I₁(k，l)+θ(5)

其中k和l是在I₁中水印嵌入位置的行数和列数，k＝22，23，24，l＝34，35，…，42，w_i＝0时θ＝0，w_i＝1时θ＝1；

步骤(7)对I′₁，I₂和I₃进行逆调和分析，利用变换公式(2)得到嵌入水印的采样信号F_t ^*；

步骤(8)对F_t ^*做反采样和反球面参数化得到加入了水印的模型M′；反采样和反球面参数化方法利用步骤(2)和(3)相应的反变换方法；

步骤(9)调整嵌入水印后的网格

通过将嵌入水印的网格进行调整来减少变形，设原始网格为M，嵌入水印后的网格为M^#，将原始网格进行球面参数化、采样、反采样和反球面参数化而不嵌入水印得到的网格为M^*，调整后的嵌入水印的网格为M^；修正公式如下：

M^＝M+α×[M^#+(M^*-M)-M]，α＜1(6)α为缩小系数；

步骤(10)提取嵌入的水印采用如下方法：

步骤a.对待检测网格M^*的X-坐标进行重新排序，使其成为由小到大上升排列方式得到一个新X-坐标的网格；

步骤b.将待检测网格M^*的X-坐标按照下式进行处理，得到网格M′：

其中，f(x)为待检测网格M^*的原始X-坐标，g(x)为处理后的待检测网格M′的原始X-坐标，T为设定的阈值，f(m)为f(x)以及前后各三个点的和

f(m)＝f(x-3)+f(x-2)+f(x-1)+f(x)+f(x+1)+f(x+2)+f(x+3)

将g(x)进行放大处理得到放大坐标f′(x)，放大函数定义如下：

f′(x)＝g(x)+(f(x)-g(x))×β    (8)

其中β为放大系数；

步骤c.利用原始网格的球面参数化的信息将网格M′的几何位置(即顶点坐标)信号F_t进行均匀采样得到采样信号F_t′，采样率为64×64，t＝1，2，3，t表示三个坐标方向；

步骤d.对F_t′做球面调和分析，得到F_t′的调和系数I_t，球面调和分析方法与嵌入水印的球面调和分析方法相同；

步骤e.将I₁中频那些嵌入水印的数值取出组成水印信息，就把水印恢复出来，即将球面调和系数中嵌入水印位置上的值取出组成序列W，将W分成四组，每组长度为6，将每组中5个最小值改为0，1个最小值改为1，则得到的编码根据水印的编码将水印信息恢复出来。