CAD工程图可逆水印方法、水印嵌入方法和水印提取方法与流程

本发明涉及一种数字水印技术领域，特别是一种CAD工程图可逆水印方法、水印嵌入方法和水印提取方法。

背景技术：

近几年，基于直方图平移的可逆水印及二维工程图可逆水印方面的研究已经取得一定的成果，主要如下：

(1)基于直方图平移的可逆水印方法

基于直方图平移的可逆水印算法一经提出便受到广泛关注。随后，G.Xuan和X.Wu相继提出了基于整数小波变换的直方图平移算法，并获得了较高的嵌入容量和水印图像质量。研究者还给出了基于预测差值的直方图平移技术，进一步改进了相关水印算法的性能。以上直方图平移算法都需要额外传输峰值点/零点信息，这在一定程度上降低了算法的实用性。现有的基于像素域的盲可逆水印方法以直方图中的最高频率像素值作为同步参考，取其邻像素值作为“峰值点”，实现直方图平移水印嵌入。其后，Tsai等人将预测差值技术与直方图平移技术结合，通过相邻像素之间的相似性构造预测误差直方图，解决了水印嵌入容量低的问题。还有采用内插误差而非预测差值的方法，通过对图像进行微小的修改而嵌入大量的水印信息。另有文献提出了一种多重预测差值和构建非对称直方图的方法，减少了水印嵌入过程中需要平移的数据个数，并结合互补嵌入策略，使两次水印嵌入前后部分数据的修改相互抵消，从而提高嵌入水印后的图像质量。

(2)二维工程图形可逆水印

对于二维工程图形的可逆水印研究，近年来也得到一定的进展。L.Li通过对线段顶点对的坐标值进行预处理操作，使其具有较强相关性后再采用差值扩大算法嵌入水印，解决了二维工程图相关性较低的问题，该算法的水印容量与载体的线段数直接相关。除此之外，现有技术提出了量化索引调制的二维CAD工程图小波域可逆水印算法，该算法首先将顶点坐标构造复值信号，然后采用改进的量化索引调制在复值信号的幅值或相角中嵌入水印信息，具有较好的不可见性，可逆性以及较高的嵌入容量，且对图形平移、旋转、均匀缩放等攻击具备较好的鲁棒性。在此基础上，研究者提出了基于小波域的二维工程图可逆水印算法，该算法先对顶点进行LPM(Log-Polar Mapping)变换和离散小波变换以获得仿射不变性，然后使用改进的量化索引调制算法将加密后的水印信息嵌入到所有的细节系数中。该算法能获得较好的不可见性和鲁棒性。现有技术还有采用改进差值扩展算法在二维CAD工程图各顶点相对距离中嵌入水印信息。改进后的算法消除了顶点间高相关性的要求，算法容易实现，且可以动态调节水印容量和不可见性。另有文献提出了两种基于预测误差直方图平移的二维CAD工程图形可逆水印算法，该算法利用二维工程图形相邻坐标值或相邻相位间的相关性，能有效的完成工程图的内容认证和隐秘通信，相比其它利用顶点相关性的可逆算法来说，在水印容量与不可见性上均具有一定优势。

因直方图平移的可逆水印方法具有对数据扰动小的优点，其在栅格图像中得到了广泛应用，然而由于二维工程图中载体数据个数少及冗余性低等特点，针对二维工程图形的直方图可逆水印的研究甚少。近年来已有很多可逆水印算法应用于二维矢量图形中，并取得一些研究成果，但是仍存在有效信息嵌入率较低，水印嵌入带来的失真较大等不足。

直方图平移的可逆水印方法通过对载体直方图平移实现高性能可逆水印嵌入。它利用直方图的零点或最小点的冗余信息，并略微修改像素的灰度值来嵌入水印信息。传统的直方图平移技术方法如下：

首先，根据载体图像生成直方图，扫描直方图中像素灰度出现频率最大的“峰值点”，记为P，找到峰值点的目的是使水印容量尽可能大，这是因为可以嵌入到图像中的比特数等于峰值点的像素个数。判断直方图中是否存在灰度频率为零的“零点”。如果存在，记零点为Z，称为有“零点情形”，如图1(a)所示。如果不存在灰度频率“零点”，称为“无零点情形”，如图1(b)。对于后者，寻找一个频率最小点并编码这些像素和像素的坐标(i,j)作为标记信息，再将其灰度值置零，则将“无零点情形”转化为“有零点情形”。

按一定顺序扫描整个图像，将位于峰值和零点之间的直方图沿零点方向平移一位，此时，原零点交换到峰值点旁，如图1(c)。再次扫描图像，遇到峰值点像素时，嵌入1bit水印。对于信息“0”，像素值保持不变；对于信息“1”，像素值变为相邻的零点(像素值加/减1)，其结果如图1(d)。

由上可知基于直方图平移可逆水印算法的水印容量取决于峰值点的个数，引入的失真一方面取决于为嵌入水印预留空间需平移的数据个数，即直方图中峰值点和零点间的元素个数；另一方面取决于嵌入信息为’1’的数目。而前者是影响图像质量的首要因素，因此在水印嵌入过程中若能减少介于峰值点和零点间数据的个数，则可有效的提高图像质量。又因为直方图平移可逆水印算法是通过对数据简单的加/减1预留空间来实现信息隐藏的，我们相信通过两次相反方向平移可使得部分的数据恢复到初值，那么图像质量也必然会有所提高。若在第一次平移过程中大部分数据集中在直方图峰值右侧，而对其左侧的数据元素向左平移(即减1)，那么在第二次直方图构建过程中，若能够使被修改的数据部分落在直方图峰值的左侧，那么此时再对其峰值右侧的数据向右平移(即加1)操作后，该数据经过两次互补修改可恢复到原始数据值。另外，由于二维CAD工程图形数据的相关性一般远低于栅格图像，如果直接采用传统的直方图平移技术，可能会导致较大的图像失真以及较低的水印容量。以往基于预测差值直方图的可逆水印大多只用一个差值来构造直方图嵌入水印信息，其直方图大致服从类拉普拉斯分布，直方图峰值两侧的差值个数基本一致，如图2。其峰值左右两侧的数据量相差为零，不具有数据分布向相反方向延伸的特点，不利于提高图像质量。

技术实现要素：

本发明所要解决的技术问题是，针对现有技术不足，提供一种CAD工程图可逆水印方法、水印嵌入方法和水印提取方法。

为解决上述技术问题，本发明所采用的技术方案是：一种CAD工程图可逆水印方法，包括水印嵌入部分、水印提取部分：

所述水印嵌入部分包括以下步骤：

1)按读取顺序提取二维工程图G中所有实体的所有顶点坐标，得到一个二维顶点坐标集合：V＝{V1,V2,...,Vi...,Vn}，其中V1,V2,...,Vi...,Vn为二维顶点坐标；遍历V中各顶点横坐标和纵坐标，获得横坐标值最小的顶点Vp和纵坐标值最小的顶点Vq为参考顶点，并将Vp和Vq在集合V中的索引位置作为提取密钥保存；其中，n为二维工程图G中所有实体的所有顶点数量；1≤p≤n；1≤q≤n；1≤i≤n，且i≠p，q；

2)将二维工程图G中所有顶点坐标转化为整数，得到二维工程图G中所有顶点相对坐标的横坐标的整数部分集合X^I和小数部分集合X^D、以及纵坐标的整数部分集合Y^I和小数部分集合Y^D或者计算Vi到参考顶点Vp和Vq间的相对相位Θi',Θi∈[0,2π]，求得Θi'的整数部分集合Θ'^I和小数部分集合Θ'^D；

3)将水印信息W＝w1,w2...wM划分为长度相等的两部分W1、W2；W1＝w1,w2...wr和W2＝wr+1,wr+2...wM；其中，1＜r＜M；M为水印信息长度；

4)对于X^I或Θ'^I中的任一元素或Θi′^I，利用公式或计算预测差值，并选取预测差值中最大者构造最大预测差值直方图E_C(max)或E_P(max)；统计E_C(max)或E_P(max)中峰值及该峰值左侧零点，分别记为Zc+、Pc+或Zp+、Pp+计算第一阶段嵌入水印后的横坐标整数值xi′^I或相对相位整数值Θi^″I：

将嵌入水印后的整数集合更新为X′^I或Θ″^I；其中，m＝n-2；

f＝1,2...r；当水印的嵌入载体是相位时，水印信息wf也相应地转化为相位大小，即wf＝1⁰或wf＝0⁰；

5)对于X′^I或Θ″^I中的任一元素xi′^I或Θi″^I，利用公式或计算预测差值，并选取预测差值中最小者构造相对坐标的最小预测差值直方图E_C(min)或相对相位的最小预测差值直方图E_P(min)，统计E_C(min)或E_P(min)中峰值及该峰值右侧零点，分别记为Zc-、Pc-或Zp-、Pp-计算第二阶段嵌入水印后的横坐标整数值xi″^I或相对相位整数值Θi″′^I：

f′＝r+1,r+2...M，当水印的嵌入载体是相位时，水印信息wf'也相应地转化为相位大小，即wf'＝1⁰或wf'＝0⁰；

6)恢复嵌入水印后的横坐标或相对相位值

其中，为X^D中的元素；P表示二维工程图G中顶点坐标小数点后的数据位数；Vpx为Vp的横坐标值；Θi′^D为Θ′^D中的元素；

7)如果待嵌入的水印比特长度大于E_C(max)和E_C(min)中峰值点的总个数，则对于Y^I中的各元素，采用与上述步骤3)～6)相同的过程进行水印嵌入；否则纵坐标值保持不变；保存嵌入水印后的顶点集合Vi^w，即得到嵌入水印后的二维CAD工程图形G^w；

所述水印提取部分包括以下步骤：

1)计算G^w中所有顶点水印嵌入后的相对坐标横坐标的整数部分集合X″^I和小数部分集合X^D；或者计算G^w中所有顶点水印嵌入后的相对相位的整数部分集合Θ″′^I和小数部分集合Θ′^D；

2)对X″^I中的任意元素xi″^I或Θ″′^I中的元素Θi″′^I，利用公式或计算最小预测差值，利用最小预测差值构造最小预测差值直方图，根据相对坐标的最小预测差值直方图峰值点pc-及零点zc-提取水印wf′′，并恢复横坐标xi′^I：

或根据相对相位的最小预测差值直方图峰值点pp-及零点zp-提取水印wf′′，并恢复相对相位整数值Θi″^I：

其中，f′＝r+1,r+2...M；

3)对元素或Θi″^I，利用公式或计算最大预测差值，根据相对坐标的最大预测差值直方图峰值点pc+及零点zc+提取水印wf′，并恢复横坐标

或根据相对坐标的最大预测差值直方图峰值点pp+及零点zp+提取水印wf′，并恢复横相对相位整数值Θi′^I：

其中，f＝1,2...r；

4)若纵坐标中也嵌入水印，对G^w中所有顶点的相对坐标集合纵坐标的整数部分采取与上述步骤1)～3)同样的方法提取水印；串接第一阶段提取的水印和第二阶段提取的水印信息，得到完整的水印信息W′＝w′1,w′2...w′M，并计算提取的水印W'与原始水印W之间的归一化相关系数NC；NC的计算公式如下：

其中，M表示水印信息W'的长度；l＝1,2,...M；

5)若NC系数为1，则根据公式或Θi＝(Θi′^I+Θi′^D)/10^P+Θq-π计算原始顶点集合V；保存恢复后的顶点集合，从而恢复出二维工程图G；其中为X^I中的元素为X^D中的元素；Θi′^D为Θ′^D中的元素。

二维工程图G中所有顶点相对坐标的计算公式如下：

其中，m＝n-2；n为二维工程图G中所有实体的所有顶点的总数量；，P≤Pmax；Pmax为二维工程图G各顶点坐标小数点后的最大数据位；Vix为Vi的横坐标值；Viy为Vi的纵坐标值；Vqy为Vq的纵坐标值；P表示水印的嵌入位置参数。该相对坐标计算公式简单，容易实现。

二维工程图G中Vi到参考顶点Vp和Vq间的相对相位Θi'的计算过程包括：计算Vp、Vq间的参考相位计算Vi到参考顶点Vp间的相位并求Θi到Θq之间的相对相位Θi'，再获取Θi'整数部分Θi'^I和小数部分Θi'^D，即得到集合Θ′^I、Θ′^D。该相对相位计算公式简单，容易实现。

本发明还提供了一种CAD工程图水印嵌入方法，包括以下步骤：

1)按读取顺序提取二维工程图G中所有实体的所有顶点坐标，得到一个二维顶点坐标集合：V＝{V1,V2,...,Vi...,Vn}，其中V1,V2,...,Vi...,Vn为二维顶点坐标；遍历V中各顶点横坐标和纵坐标，获得横坐标值最小的顶点Vp和纵坐标值最小的顶点Vq为参考顶点，并将Vp和Vq在集合V中的索引位置作为提取密钥保存；其中，n为二维工程图G中所有实体的所有顶点数量；1≤p≤n；1≤q≤n；1≤i≤n，且i≠p，q；

2)将二维工程图G中所有顶点坐标转化为整数，得到二维工程图G中所有顶点相对坐标的横坐标的整数部分集合X^I和小数部分集合X^D、以及纵坐标的整数部分集合Y^I和小数部分集合Y^D或者计算Vi到参考顶点Vp和Vq间的相对相位Θi',Θi′∈[0,2π]，求得Θi'的整数部分集合Θ'^I和小数部分集合Θ'^D；

3)将水印信息W＝w1,w2...wM划分为长度相等的两部分W1、W2；W1＝w1,w2...wr和W2＝wr+1,wr+2...wM；其中，1＜r＜M；M为水印信息长度；

4)对于X^I或Θ'^I中的任一元素或Θi′^I，利用公式或计算预测差值，并选取预测差值中最大者构造最大预测差值直方图E_C(max)或E_P(max)；统计E_C(max)或E_P(max)中峰值及该峰值左侧零点，分别记为Zc+、Pc+或Zp+、Pp+计算第一阶段嵌入水印后的横坐标整数值xi′^I或相对相位整数值Θi″^I：

将嵌入水印后的整数集合更新为X′^I或Θ″^I；其中，m＝n-2；当水印的嵌入载体是相位时，水印信息wf也相应地转化为相位大小，即wf＝1⁰或wf＝0⁰；

5)对于X′^I或Θ″^I中的任一元素xi′^I或Θi″^I，利用公式或计算预测差值，并选取预测差值中最小者构造相对坐标的最小预测差值直方图E_C(min)或相对相位的最小预测差值直方图E_P(min)，统计E_C(min)或E_P(min)中峰值及该峰值右侧零点，分别记为Zc-、Pc-或Zp-、Pp-计算第二阶段嵌入水印后的横坐标整数值xi^″I或相对相位整数值Θi″′^I：

当水印的嵌入载体是相位时，水印信息wf'也相应地转化为相位大小，即wf'＝1⁰或wf'＝0⁰；

6)恢复嵌入水印后的横坐标或相对相位值

其中，为X^D中的元素；P表示二维工程图G中顶点坐标小数点后的数据位数；Vpx为Vp的横坐标值；Θi′^D为Θ′^D中的元素；

7)如果待嵌入的水印比特长度大于E_C(max)和E_C(min)中峰值点的总个数，则对于Y^I中的各元素，采用与上述步骤3)～6)相同的过程进行水印嵌入；否则纵坐标值保持不变；保存嵌入水印后的顶点集合Vi^w，即得到嵌入水印后的二维CAD工程图形G^w。

相应地，本发明的CAD工程图水印提取方法包括以下步骤：

1)计算上述水印嵌入后得到的二维CAD工程图形G^w中所有顶点水印嵌入后的相对坐标横坐标的整数部分集合X″^I和小数部分集合X^D；或者计算G^w中所有顶点水印嵌入后的相对相位的整数部分集合Θ″′^I和小数部分集合Θ′^D；

2)对X″^I中的任意元素xi″^I或Θ″′^I中的元素Θi″′^I，用公式或计算最小预测差值，利用最小预测差值构造最小预测差值直方图，根据相对坐标的最小预测差值直方图峰值点pc-及零点zc-提取水印wf′′，并恢复横坐标xi′^I：

或根据相对坐标的最小预测差值直方图峰值点pp-及零点zp-提取水印wf′′，并恢复相对相位整数值Θi″^I：

3)对元素或Θi″^I，利用公式或计算最大预测差值，根据相对坐标的最大预测差值直方图峰值点pc+及零点zc+提取水印wf′，并恢复横坐标

或根据相对相位的最大预测差值直方图峰值点pp+及零点zp+提取水印wf′，并恢复相对相位整数值Θi′^I：

4)若纵坐标中也嵌入水印，对G^w中所有顶点的相对坐标集合纵坐标的整数部分Y″^I采取与上述步骤1)～3)同样的方法提取水印；串接第一阶段提取的水印和第二阶段提取的水印信息，得到完整的水印信息W′＝w′1，w′2...w′M，并计算提取的水印W'与原始水印W之间的归一化相关系数NC；NC的计算公式如下：

其中，M表示水印信息W'的长度；l＝1,2,...M；

5)若NC系数为1，则根据公式或Θi＝(Θi′^I+Θi′^D)/10^P+Θq-π计算原始顶点集合V；保存恢复后的顶点集合，从而恢复出二维工程图G；其中为X^I中的元素；为X^D中的元素；Θi′^D为Θ′^D中的元素。

与现有技术相比，本发明所具有的有益效果为：本发明对工程图形的当前顶点的整数部分采用多重预测的方法，计算多个预测差值，并通过选取合适的差值分别构造非对称直方图，使得峰值左右两侧的数据量相差悬殊且具有数据分布向相反方向延伸的特点。在水印嵌入过程中，对两个直方图所对应的顶点数据进行反向修改，使得两次水印嵌入部分的顶点数据恢复到初始值，从而提高图像质量。对二维工程图而言，其部分相邻顶点对或周边的顶点具有相同的横坐标或纵坐标，或者具有相同的坐标值差值，也就是说，通过周边的几个预测差值都能准确地预测目标顶点，因而本发明方法能够有效降低水印嵌入过程中需平移的数据个数，在水印容量与不可见性上具有一定优势，能有效的完成工程图的内容认证和隐秘通信；具备严格的可逆性，在图形未经恶意攻击时，原始的二维工程图能够在水印提取后完全恢复。

附图说明

图1(a)为“零点情形”示意图；图1(b)为“无零点情形”示意图；图1(c)为原零点交换到峰值点旁示意图；图1(d)为信息“1”像素值变为相邻的零点示意图。

图2为传统预测差值直方图；

图3(a)为本发明Emax直方图；图3(b)为本发明Emin直方图；

图4(a)和图4(b)所示为非对称预测差值直方图平移水印嵌入和提取水印的流程图；

图5(a)～图5(i)为本发明原始图形、嵌入水印图形及恢复图形与对应预测误差直方图；

图6为Asym_HSC算法选取不同的嵌入位置与嵌入容量的关系图；

图7(a)为P＝3的预测差值直方图；图7(b)为P＝4的预测差值直方图；图7(c)为P＝5的预测差值直方图；

图8为Asym_HSC不可见性与嵌入位置关系图。

具体实施方式

本发明方法具体包括：

Step1:构造非对称直方图

首先，对待处理数据序列按(1)、(2)、(3)式计算预测差值：

当i＝1时，其预测差值可计算为：

eci1＝eci2＝xi. (1)

当i＝2时，其预测差值可计算为：

当i>＝3时，其预测差值可计算为：:

选取预测差值中较大者eci(max)＝max(eci1,eci2)构造最大预测差值直方图EC(max)＝ec1(max),...ecm(max)如图3(a)，不同于一般的预测误差直方图服从类拉普拉斯分布，本发明通过选取两个预测差值中较大者作为构造直方图的元素，则大部分差值数据集中在峰值点的右侧，而峰值左侧的元素个数要远少于其右侧元素个数，此时直方图呈不对称分布趋势。

同理，可选取较小的预测差值eci(min)＝min(eci1,eci2)构造直方图EC(min)＝ec1(min),...ecm(min)，如图3(b)所示，此时直方图大部分数据集中在峰值点左侧。

Step2:双重互补水印嵌入

以X坐标为例，先构造最大预测差值直方图Emax后构造最小预测差值直方图Emin阐述本发明的方法。，在Emax直方图中确定峰值Pc+及其左侧零点Zc+，按顺序扫描所有的预测差值，水印嵌入过程为：当预测差值在区间(Zc+,Pc+)内时，该预测差值对应的顶点数据向左移动一位为水印嵌入预留空间。遇到预测差值峰值点时，其对应的顶点数据整数部分嵌入1bit水印：对于信息“0”，其值保持不变；对于信息“1”，其值减1。以X坐标为例，其水印嵌入过程如下：

确定Emin中峰值P-及右侧零点Z-，再次按顺序扫描所有的预测差值，若其在区间(P-，Zc-)范围内，该预测差值对应的顶点整数部分向右移动一位为水印嵌入预留空间。若预测差值正好为峰值点时，其对应的值嵌入1bit水印：对于信息“0”，保持不变；对于信息“1”，其值加1。其水印嵌入过程如下：

由于直方图的峰值和零点在本发明方法中具有重要位置，下面将对直方图中峰值和零点的选择作详细阐述。

定义1：直方图零点。在预测差值直方图中，如果差值ei和其所发生的次数(频数)f(ei)存在这样的关系：若f(ei)＝0，称此时的ei为直方图上的一个零点。

定义2：直方图峰值。在预测差值直方图中，如果差值ei和其所发生的次数(频数)f(ei)存在这样的关系：若f(ei)＝max{f(e1),f(e2),...f(em)}，则称ei为差值直方图上的峰值。

本发明利用非对称预测差值直方图双重互补嵌入的方法分别对坐标预测差值和相位预测差值进行水印嵌入，为方便描述，定义在相对坐标差值中嵌入水印的算法为Asym_HSC(本发明实施例以Asym_HSC为例进行说明)，在相对相位差值中嵌入水印的算法为Asym_HSP。通过对应的预测差值修改顶点数据和相位数据，避免了因实际值与预测值差异较大而引入较大失真。图4和图5所示为非对称预测差值直方图平移水印嵌入和提取水印的流程图。对于嵌入部分，水印信息W是二进制比特流，水印嵌入分互补的两个阶段，分别为最大差值直方图水印嵌入和最小差值直方图水印嵌入。预测差值直方图的创建是为修改顶点或相位数据及嵌入水印提供依据。对于水印提取部分，采用与水印嵌入部分相同的变换则可以提取水印并恢复原始图形。主要可以分为三个阶段：数据预处理、水印嵌入和水印提取。

1.数据预处理

这里用一维顶点向量V表示从所有实体中提取的二维顶点坐标数据的集合，其中：V＝{(V1x,V1y),(V2x,V2y),...,(Vix,Viy),...,(Vnx,Vny)}，其中(Vix,Viy)为第i个顶点对应的横坐标和纵坐标，n为顶点的总数量。

(1)相对坐标计算

Step1：遍历二维CAD工程图各顶点坐标值，分别选取横坐标值和纵坐标值最小的顶点作为参考顶点Vp和Vq。

Step2：计算相对坐标序列：通常二维CAD工程图的坐标由固定精度(6～8位)的浮点数组成。为了适用于直方图平移技术，所有的顶点坐标需要先转化为整数。假设小数点后的最大数据位为Pmax，则整数坐标的计算公式如式(6)所示。

其中P表示取小数点后数据位数，floor(·)表示取整数部分。依此可计算得到一个相对坐标的整数集合X^I,Y^I和小数集合X^D,Y^D。

其中m＝n-2(当p＝q时，则m＝n-1)。X^I,Y^I为Asym_HSC差值直方图统计的载体数据，可以通过直方图平移和修改嵌入水印信息。

(2)相对相位计算

Step1：遍历二维CAD工程图各顶点横坐标值，选取横坐标值最小的顶点作为参考顶点Vp。并使用密钥K随机选取工程图顶点集中另一顶点Vq作为第二个参考顶点。

Step2：计算相对相位序列：首先计算两个参考顶点间的相位对V中的任意顶点Vi(i≠p，q)，计算Vi与Vp之间幅值和相位于是，可计算得到Vi到参考顶点Vp和Vq间的相对相位：

Θi'＝Θi-Θq+π，Θi′∈[0,2π]. (8)

依此由顶点集合V可计算得到一个相对相位集合Θ′＝{Θ1'，Θ2′，...Θm'}，其中m＝n-2。该相位集合Θ′即为Asym_HSP算法的载体数据。

Step3：对于相对相位集Θ'中各个元素，再采用式(9)求整数部分Θi'^I和小数部分Θi'^D。

相应的逆变换为：

Θi＝(Θi'^I+Θi'^D)/10^P+Θq-π； (10)

2水印嵌入

Step1：按读取顺序提取二维工程图G中所有实体的顶点数据，得到一个二维顶点集合：V＝{V1,V2,...,Vi...,Vn}。遍历V中各顶点横坐标和纵坐标，获得参考顶点Vp和Vq，并将其在集合V中的索引位置作为提取密钥保存。

Step2：根据上述数据预处理计算方法，计算整数集X^I和小数集X^D。

Step3：为适应双重嵌入，将水印信息W划为W1和W2两部分，在水印嵌入的第一阶段将W1隐藏到顶点数据中，水印嵌入第二阶段将W2隐藏到顶点数据中。

Step4：对于整数集X^I中各个元素，先采用式(1)计算两个邻近差值，选取较大者构造最大差值直方图，并统计差值直方图中峰值和及其左侧零点，分别记为Z+和P+(一般峰值为ei＝0处)。对预测差值在Z+与P+之间的数据值向左平移一个单位，其数据修改和水印嵌入如下：

其中xi′^I为第一阶段嵌入水印后的坐标整数值。

Step5：整理第一阶段嵌入水印后的顶点数据，继续采用(3)计算邻近差值，选取较小者构造最小差值直方图，并统计差值直方图中峰值和及其右侧零点，分别记为Zc-和Pc-(一般峰值为ei＝0处)。对预测差值在Zc-与Pc-之间的数据值向右平移一个单位，其数据修改和水印嵌入如下：

其中xi″^I为第二阶段嵌入水印后的坐标整数值。

Step6：恢复嵌入水印后的横坐标：

Step7：如果待嵌入的水印比特长度大于Emax直方图和Emin直方图中峰值点的总个数，则对于Y^I中的各元素，采用与X^I相同的方式处理，否则纵坐标值保持不变。保存嵌入水印顶点集合Vi^w，即可得到嵌入水印后的二维CAD工程图形G^w。

Pmax表示工程图坐标点数值小数点后的最大有效位，P表示水印的嵌入位置，即在第P个小数点后嵌入。假设工程图的一个坐标点大小为6253.123456，则Pmax＝6，若取P＝1，即在小数点后的第一位嵌入。

3水印提取

为正确提取水印信息，提取部分的操作和嵌入部分的操作顺序相反，若嵌入过程中先构造E(max)后构造E(min)嵌入水印信息，那么在提取水印时应先构造E(min)提取水印再构造E(max)提取水印和恢复原始数据，反之亦然。具体操作如下：

Step1：根据嵌入过程的索引值获取参考顶点Vp，再计算含水印工程图G^w中所有顶点的相对坐标集合的整数部分和小数部分X′^I和X′^D。

Step2：对X″^I采用(3)进行邻近最小差值预测，逐步提取水印和恢复数据如式(14)、(15)，完成第一阶段水印的提取。

Step3：继续对X^I'坐标采用(2)进行邻近最大差值预测，提取水印和恢复数据如(16)、(17)式，完成第二阶段的水印提取。

Step4：若Y坐标中也嵌入水印，对Y^I”可采取与X^I”同样的方法提取水印和恢复数据。计算提取水印信息W'与原始水印W之间的归一化相关系数NC。相关系数的计算公式如下：

Step5：若NC系数为1.0，根据公式(12)，计算得到原始顶点集合V。保存恢复后的顶点序列，从而恢复出原始二维CAD工程图G。

至此，水印及工程图数据已完全提取恢复和恢复，同样的方法可提取二维CAD工程图相对相位中的水印信息。

实验结果与分析

为验证和分析本发明方法的有效性，在配置为CPU 3.0GHz，RAM 2.0GB，Windows XP，Visual C++6.0，DWGDirect 4.0，AutoCAD 2008，以及Matlab R2010a的个人计算机上对50幅标准二维CAD工程图形做仿真实验，这里采用矢量图中常用指标RMSE(均方根误差)来量化衡量原始图形与水印图形之间的差异，RMSE的计算公式如下：

其中N为顶点总数，V和V^w分别为原始顶点与嵌入水印后顶点。

图5(a)～图5(i)为工程图G2采用Asym_HSC方法在P＝1时嵌入水印和可逆恢复后的图形及对应预测误差直方图，其中5(d)为选取预测差中最大差值Emax构造的直方图，5(g)为选取预测差中最小差值Emin构造的直方图，可以看出其预测差值直方图不再服从类拉普拉斯分布，而是呈非对称分布，在Emax直方图中，峰值左侧的数据明显少于其右侧的数据，并且该直方图拥有明显的峰值，嵌入水印时对其左侧预测值对应的顶点数据进行修改，在保证较大嵌入容量的同时减少图像失真。5(e)、5(h)均为经数据平移和修改后的误差直方图，由图可知，其直方图形状基本保持不变，但峰值存在一定程度的下降，这主要是由于水印嵌入过程导致峰值处的差值向相邻差值处转移。恢复后图形的预测误差直方图5(f)、5(i)，且具有原始图形完全一致的直方图形状，由此也表明G2已无损恢复。

通常，只要恢复的二维CAD工程图顶点坐标在小数点后Pmax位有效数字与原始图形完全一致，则表示二者是相同的，即水印算法是可无损恢复的。50幅二维CAD工程图形的基本信息以及使用Asym_HSC和Asym_HSP算法满嵌水印，完全提取水印后恢复的图形与原始图形的RMSE值如表1所示。

表1原始图形与恢复图形的RMSM

50幅工程图中恢复图形与原始图形之间的RMSE值范围在10^-13～10^-16之间，由此表明本文算法具有良好的可逆性。此外，由图5(c)可见，可逆恢复后的工程图形与原始图形从视觉上完全一致。

在水印信息的嵌入过程中，只有当扫描图形遇到各直方图的峰值点即Emin(X)峰值P(x)-、Emax(X)峰值P(x)+和Emin(Y)峰值P(y)-、Emax(Y)峰值P(y)+才会嵌入1bit水印信息，故可由直方图的P(x)-和P(x)+、P(y)-、P(y)+精确计算嵌入容量Capacity：

Capacity＝(f(P(x)-)+f(P(x)+)f(P(y)-)+f(P(y)+))/N (20)

这里采用bpv(bit per vertex)用来度量。

针对文献(1、D.Wu,G.Z.Wang.Reversible Watermarking of SVG Graphics.In:Proceedings of 2009International Conference on Communications and Mobile Computing,2009:385-390；2、L.Li,Q.Li,X.Fang,B.Qian.Reversible watermarking algorithm based on difference expansion for 2D engineering drawings.Journal of Image and Graphics,2010,15(2):372-376；3、F.Peng,Liu Y,Long M.Reversible watermarking for 2D CAD engineering graphics based on improved histogram shifting[J].Computer-Aided Design,2014,49:42-50)提出的可逆水印方法也测试了其水印容量，如表2所示，水印容量高于其他方法。由表中数据可见，当设置水印嵌入位置相同时，Asym_HSC和Asym_HSP算法的水印容量均优于文献(F.Peng,Liu Y,Long M.Reversible watermarking for 2D CAD engineering graphics based on improved histogram shifting[J].Computer-Aided Design,2014,49:42-50)的同类方法，

表2几种不同方法水印容量比(单位：比特/顶点)

图6所示为Asym_HSC算法选取不同的嵌入位置与嵌入容量的关系：随着嵌入位置的后移，水印容量有一定程度的降低，之后逐渐趋于稳定。以工程图G2为例，其在不同嵌入位置下的差值直方图如图7(a)、图7(b)、图7(c)所示。当P＝3,4,5时的差值直方图均保持相对稳定的形状，且峰值也基本保持不变，由此表明当P>3时水印容量的稳定性。

本发明对不同可逆水印方案(1、D.Wu,G.Z.Wang.Reversible Watermarking of SVG Graphics.In:Proceedings of 2009International Conference on Communications and Mobile Computing,2009:385-390；2、L.Li,Q.Li,X.Fang,B.Qian.Reversible watermarking algorithm based on difference expansion for 2D engineering drawings.Journal of Image and Graphics,2010,15(2):372-376；3、F.Peng,Liu Y,Long M.Reversible watermarking for 2D CAD engineering graphics based on improved histogram shifting[J].Computer-Aided Design,2014,49:42-50)的不可见性进行了对比分析，另外，由于文献3与本实验都是基于直方图的可逆水印算法，为实现算法的可比性，本发明设置p＝1和嵌入容量相同的情况下中对这两种的方法进行了比较，实验结果如表3所示。

表3原始图形与水印图形的RMSE值比较

由表3可以看出，Asym_HSC算法的不可见性要优于同类算法，由于Asym_HSC算法直接对横坐标值和纵坐标值处理，在最极端的情况下(即横坐标和纵坐标同时平移一个单位)，坐标值的移动对该顶点的变化量至多为

因此具有较好的不可见性；而Asym_HSP算法中对相位值10^-P的改变量，而对顶点值的改变值相应为：

由式(20)、(21)可知，当顶点坐标的幅值大于时，通过坐标嵌入水印比通过相位嵌入水印引入的失真小。表4是对50幅工程图形顶点中幅值Λi大于及小于等于分布的统计，从表中可以看出，相对于Asym_HSC算法而言，Asym_HSP算法对图形造成的扰动更大。就同类算法而言，Asmy_HSC和Asym_HSP与文献3的HS_C和HS_P相比均具有一定的优势，其原因是通过构造非对称直方图，减少了需平移的顶点数目，从而提高水印后图像的质量。并且在互补嵌入水印过程中，部分顶点经过两次相反方向的平移后误差相互抵消，可进一步提高水印不可见性。

表4工程图形顶点幅值对比

图8所示为Asmy_HSC算法不可见性与嵌入位置的关系图。由该图可分析：随着嵌入位置的后移，RMSE值显著降低。因此为了确保水印的不可见性和较大的水印容量，通常可以将嵌入位置参数设置为1～3。

为了验证本发明方法的鲁棒性能，对上述4幅标准工程图形在采用算法Asym_HSC和Asym_HSP嵌入水印后，分别施加平移、旋转、均匀缩放操作，并计算提取的水印信息与原始水印信息的相关系数，结果如表5所示。Asym_HSC算法对于平移操作的鲁棒性较好，而对旋转和缩放操作的则相当脆弱；而Asym_HSP算法对于这三种操作的NC均为1.0，即Asym_HSP算法能够抵抗平移、旋转和均匀缩放操作。

表5常规操作后提取的水印信息与原始水印的相关系数