圖像數(shù)字水印+matlab程序

1、第三章 圖像數(shù)字水印的方案3.1 圖像數(shù)字水印的技術(shù)方案在數(shù)據(jù)庫中存儲在國際互聯(lián)網(wǎng)上傳輸?shù)乃D像一般會被壓縮，有時達(dá)到很高的壓縮比。因此，數(shù)字水印算法所面臨的第一個考驗(yàn)就是壓縮。JPEG和EZW（Embedded Zero-Tree Wavelet）壓縮是最常見的兩種壓縮方法。JPEG是基于離散余弦變換域的壓縮方法，而EZW是基于小波變換域的壓縮方法。前人的研究證明采用與壓縮算法相同的變換域水印方法，對于壓縮的穩(wěn)健性較強(qiáng)。因此，我研究圖像文件水印算法主要集中在變換域算法及利用人眼視覺特性上。數(shù)字水印的嵌入要求即要考慮視覺透明性，又要保證嵌入水印后圖像的穩(wěn)健性，這兩個方面存在著矛盾。保證視覺透

2、明性，就要將水印嵌入到人眼不敏感區(qū)，也就是嵌入到圖像的高頻分量中。而多數(shù)圖像處理方法對于圖像高頻部分的損壞程度較高，如有損壓縮、高頻濾波等。水印很容易在經(jīng)歷圖像處理的過程中丟失。這樣，則無法保證圖像數(shù)字水印的穩(wěn)健性。如果要獲得很好的穩(wěn)健性，數(shù)字水印應(yīng)加在人眼敏感的低頻部分，圖像的大部分能量集中在低頻部分，如果對于低頻部分進(jìn)行處理，水印固然會失去，而圖像也沒有了利用價值，然而，水印的嵌入會對圖像的質(zhì)量有非常大的影響，這又無法保證視覺透明性。數(shù)字水印算法的實(shí)現(xiàn)基本分為三個部分：宿主圖像的變換，水印的嵌入和水印的檢測，分別描述如下。3.2 基于DCT域的圖像數(shù)字水印技術(shù)離散余弦變換（Discrete

3、 Cosine Transform）屬于正交變換圖像編碼方法中的一種。正交變換圖像編碼始于1968年。當(dāng)時安德魯斯（Andrews）等人發(fā)現(xiàn)大多數(shù)自然圖像的高頻分量相對幅度較低，可完全舍棄或者只用少數(shù)碼字編碼，提出不對圖像本身編碼，只對其二維傅立葉（DFT）系數(shù)進(jìn)行編碼和傳輸。但DFT是一種正交變換，運(yùn)算量很大，常常使實(shí)時處理發(fā)生困難，第二年他們就用Walsh-Hadamard變換（WHT）取代DFT可以使運(yùn)算量明顯減少，這是因?yàn)閃HT變換只有加減法而無需乘法。但是更有意義的是離散余弦變換和離散正旋變換的出現(xiàn)，它們具有快速算法，精確度高。其中最重要的是1974年提出的DCT，因?yàn)槠渥儞Q矩陣的基

4、向量很近似于托伯利茲矩陣的特征向量，而托伯利茲矩陣又體現(xiàn)了人類語言及圖像信號的相關(guān)性。因此，DCT常常被認(rèn)為是語音與圖像信號變換的準(zhǔn)最佳變換。圖像是二維的，所以在研究時主要用到二維DCT，以及二維IDCT來對圖像進(jìn)行處理。3.2.1 離散余弦變換（DCT）的定義數(shù)字圖像X（m,n）是具有M行N列的一個矩陣。為了同時減弱或去除圖像數(shù)據(jù)相關(guān)性，可以運(yùn)用二維DCT，將圖像從空間域轉(zhuǎn)換到DCT變換域。根據(jù)定義，二維離散余弦變換（DCT）定義如下：式中：m，k=0,1,M-1; n， =0,1,N-1。 二維逆離散余弦變換（IDCT）的定義如下：式中：m，k=0,1,M-1; n，=0,1,N-1。3.

5、2.2 離散余弦變換的特點(diǎn)在基于DCT的變換編碼中，圖像是先經(jīng)分塊（8×8或16×16）后再經(jīng)DCT，這種變換是局部的，只反映了圖像某一部分的信息。當(dāng)然也可以對整幅圖像的特點(diǎn)，但是運(yùn)算速度比分塊DCT要慢。圖像經(jīng)DCT后，得到的DCT圖像有三個特點(diǎn)：一是系數(shù)值全部集中到0值附近（從直方圖統(tǒng)計(jì)的意義上），動態(tài)范圍很小，這說明用較小的量化比特?cái)?shù)即可表示DCT系數(shù)；二是DCT變換后圖像能量集中在圖像的低頻部分，即DCT圖像中不為零的系數(shù)大部分集中在一起（左上角），因此編碼效率很高。三是沒有保留原圖像塊的精細(xì)結(jié)構(gòu)，從中反映不了原圖像塊的邊緣、輪廓等信息，這一特點(diǎn)是由DCT缺乏時局域

6、性造成的。如下左圖31是原始圖像經(jīng)過DCT變換后的系數(shù)圖像為圖32。兩條線劃分出圖像的低頻、中頻和高頻分別所在的矩形區(qū)域?？梢钥闯?，圖像DCT變換后大部分參數(shù)接近于零，只有左上角的低頻部分有較大的數(shù)值，中頻部分參數(shù)值相對較小，而大部分高頻參數(shù)值非常小，接近于零。. 圖 31：原圖像圖 32：變換后的系數(shù)圖像3.2.3 離散余弦變換的數(shù)字水印算法根據(jù)離散余弦變換后的參數(shù)性質(zhì)，本文采用了以ZigZag方式重排變換域系數(shù)的方法，選出中頻分量，用數(shù)字水印序列對其進(jìn)行非線性調(diào)制。水印檢測時，待檢測圖像仍按比方式選擇變換域系數(shù)，與待水印進(jìn)行相關(guān)運(yùn)算，與閾值比較來判斷是否所含水印。離散余弦域的數(shù)字水印算法的

7、具體實(shí)現(xiàn)分為三步：宿主圖像的變換，數(shù)字水印的嵌入，數(shù)字水印的檢測。3.2.3.1 宿主圖像的DCT變換對于N×N大小的256灰度級的宿主圖像I進(jìn)行N×N二維離散余弦變換（DCT）。以ZigZag方式對于DCT變換后的圖像頻率系數(shù)重新排列成一維向量Y=y1, y2,yN×N.并取出序列中第L+1到L+M的中頻系數(shù)部分,得到Y(jié)L= YL+1, YL+2, YL+M 3.2.3.2 數(shù)字水印的嵌入假設(shè)數(shù)字水印W為一服從標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)實(shí)數(shù)序列，用數(shù)字序列表示為W=W1 ，W2 ，WM 。用W對Y序列中第L+1到L+M的中頻系數(shù)部分的值進(jìn)行修改，按以下公式進(jìn)行： 經(jīng)過修

8、改的系數(shù)序列Y = Y1, Y2,. YN×N 以ZigZag逆變換形式重組，再進(jìn)行N×N DCT逆變換，得到嵌有數(shù)字水印的圖像I。3.2.3.3 數(shù)字水印的檢測待檢測的可能含有水印的圖像I" 。假設(shè)I"未損失大量信息，可以近似認(rèn)為I"= I。在此假設(shè)下可以運(yùn)用統(tǒng)計(jì)的方法來檢測水印。（1）待檢水印域待檢圖像中頻系數(shù)相關(guān)性的測定同樣對I進(jìn)行DCT變換，以ZigZag方式將DCT系數(shù)排成一維向量Y "= Y1 ", Y2 ",. YN×N"。由于假設(shè)I"=I，則Y"= Y 。取出

9、Y"（等于Y）中第L+1到L+M的中頻系數(shù)部分YL "= YL+1 ", YL+2",. YL+M"。假設(shè)待檢測的數(shù)字水印X=X1，X2，. XM為一符合標(biāo)準(zhǔn)正態(tài)分布的實(shí)數(shù)偽隨機(jī)序列。則可以通過待檢水印與圖像中頻系數(shù)作相關(guān)運(yùn)算來判斷是否所加入了水印。只有在待檢水印為所加入的水印時，才能得到較大的相關(guān)值。否則相關(guān)值很小，接近于零。 用符號E表示數(shù)學(xué)期望，得到： （2）閾值的確定根據(jù)中心極限定理，參照水印匹配與不匹配兩種情況得到閾值為。由于原始圖像難以得到，因此從實(shí)用性出發(fā)，閾值定義為： 綜上所述，滿足時，則表明檢測到匹配水印。否則，未檢測到匹配的

10、水印。3.3 MATLAB工具簡介3.3.1. 簡介Matlab是當(dāng)前在國內(nèi)外十分流行的工程設(shè)計(jì)和系統(tǒng)仿真軟件包。它是MathWorks公司于1982年推出的一套高性能的數(shù)值計(jì)算和可視化軟件，它集數(shù)值分析、矩陣運(yùn)算、信號處理和圖形顯示于一體，構(gòu)成了一人方便的、界面友好的用戶環(huán)境。Matlab的推出得到了各個領(lǐng)域?qū)＜摇W(xué)者的廣泛關(guān)注，其強(qiáng)大的擴(kuò)展功能為各個領(lǐng)域的應(yīng)用提供了基礎(chǔ)。由各個專家學(xué)者相繼推出了MATLAB工具箱，其中的信號處理(signal processing)、控制系統(tǒng)(control system)、神經(jīng)網(wǎng)絡(luò)(neural network)、圖像處理(image processin

11、g)、魯棒控制(robust control)、非線性系統(tǒng)控制設(shè)計(jì)(nonlinear system control design)、系統(tǒng)辨識(system identification)、最優(yōu)化(optimization)、模糊邏輯(fuzzy logic)、小波(wavelet)、通信(communication)、統(tǒng)計(jì)(statistics)等工具箱，這些工具箱給各個領(lǐng)域的研究和工程應(yīng)用提供了有力的工具，借助于這些“巨人肩上的工具”，各個層次的研究人員可直觀、方便地進(jìn)行分析、計(jì)算及設(shè)計(jì)工作，從而大大地節(jié)省了時間。3.3.2. MATLAB研究數(shù)字水印的優(yōu)點(diǎn) 集成了DCT、DWT等函數(shù)有豐

12、富的小波函數(shù)和處理函數(shù)，這不僅方便了研究人員，而且使源程序簡潔明了、易實(shí)現(xiàn)。 強(qiáng)大的數(shù)學(xué)運(yùn)算功能。能夠方便、高效地實(shí)現(xiàn)音頻、視頻中的大量矩陣運(yùn)算。 提供了圖像處理工具箱、小波分析工具箱、數(shù)字信號處理工具箱。用來編制跨數(shù)字圖像處理技術(shù)、數(shù)字信號處理等多學(xué)科的數(shù)字水印技術(shù)是非常好的選擇。 MATLAB與目前最強(qiáng)大的編程工具Visual C+具有良好的接口。3.3.3. MATLAB函數(shù)介紹在介紹函數(shù)之前，我們必須明確一點(diǎn)：作水印程序時，處理的圖像數(shù)據(jù)是二維信號，而聲音信號是一維信號。這里，我們僅僅簡單介紹與水印有關(guān)的函數(shù)。 數(shù)據(jù)輸入輸出函數(shù)imread（）和imwrite（）：可以讀寫bmp,jp

13、g/jpeg, tif/tiff, png, hdf, pcx, wxd格式文件。讀索引文件時，還可以得到相應(yīng)的調(diào)色板數(shù)據(jù)。auread（）、auwrite（）、wavread（）和wavwrite（）：可以方便地讀寫au和wav文件，并可控制其中的位及頻率。 圖像顯示imshow()：顯示一幅圖像；imfinfo()：可以得到讀入圖像的信息。如文件的大小、格式、格式版本號、圖像的高度、寬度、顏色類型（真彩色，灰度圖還是索引圖）等。 變換頻函數(shù)對信號采用不同的變換，是實(shí)現(xiàn)頻域法水印的至關(guān)重要的一步，MATLAB中對一維信號和二維信號分別提供了各種變換和逆變換函數(shù)。u 離散余弦變換（DCT）dc

14、t(),dct2():分別實(shí)現(xiàn)一維信號和二維信號的DCT（離散余弦變換）；idct(),idct2():分別實(shí)現(xiàn)一維信號和二維信號的IDCT（逆向離散余弦變換）；u 離散小波變換（DWT）dwt(),dwt2():分別實(shí)現(xiàn)一維信號和二維信號的DWT（離散小波變換）idwt(),idwt2():分別實(shí)現(xiàn)一維信號和二維信號的IDWT（離散小波變換）；Wavedec2():多級二維小波分解函數(shù)；Waveinfo():提供小波包中所有的小波信息；攻擊函數(shù)對算法進(jìn)行攻擊測試是對水印魯棒性檢測的一種重要手段，一個好的水印算法必須經(jīng)過各種攻擊測試才能對之做出客觀的評價。MATLAB中的許多函數(shù)可以直接用來做

15、攻擊測試。旋轉(zhuǎn)：rotate（）可以對圖像進(jìn)行任意角度的旋轉(zhuǎn)；剪裁：imcrop()可以按精確定位的各點(diǎn)坐標(biāo)進(jìn)行剪裁；濾波：filter()和filter2()可實(shí)現(xiàn)對一維信號和二維信號的濾波；抖動：dither()對圖像進(jìn)行抖動；抖動攻擊考驗(yàn)水印魯棒性的一個很好的攻擊；jpeg壓縮：imwrite()中jpg和quality參數(shù)能對圖像進(jìn)行可控jpg壓縮；加各種噪聲：imnoise()可以對圖像加入各種噪聲，如白噪聲、椒鹽噪聲等，加入噪聲是對水印魯棒性考驗(yàn)的一種常見的攻擊；放大/縮?。篿mresize()可以以指定的插值方法來對圖像進(jìn)行放大和縮小。第四章： 圖像數(shù)字水印技術(shù)的實(shí)現(xiàn)4.1 基于

16、離散余弦變法（DCT）實(shí)現(xiàn)數(shù)字水印技術(shù) 打開原始及水印圖像：subplot(2,2,1)I=uigetfile('*.bmp','打開原始彩色圖像文件');RGB=imread(I);image(RGB);title('原始彩色圖像');subplot(2,2,2)I=uigetfile('*.bmp','打開水印灰度圖像文件');imshow(I);title('灰度水印圖像');subplot(2,2,3)H=imread(I);J=dct2(H);imshow(log(abs(J),),col

17、orbar;title('水印圖像經(jīng)DCT變換后能量分布情況')運(yùn)行結(jié)果： 水印全過程：0%水印加入程序Q=input('請輸入放縮因子的值(建議小于1)：Q')subplot(2,3,1)RGB=imread('浙江臺州學(xué)院','jpg');imshow(RGB);title('原始圖像');subplot(2,3,2)N=dct2(RGB(:,:,3);imshow(log(abs(N),),colorbar;title('Y分量能量分布');subplot(2,3,4)I=imread(

18、9;lena1','bmp');imshow(I);title('灰度水印圖像');subplot(2,3,5)M=dct2(I);imshow(log(abs(M),),colorbar;title('水印能量分布');subplot(2,3,6)J=M(1:128,1:128);J(128:364,128:400)=0;J=rot90(J);J=rot90(J);J(365:600,401:750)=0;J=rot90(J);J=rot90(J);N=N+Q*J;K=idct2(N);RGB(:,:,3)=K;imshow(RGB);

19、title('加入水印后圖像');%水印提取程序subplot(2,3,3)RGB1=imread('浙江臺州學(xué)院','jpg');N=dct2(RGB(:,:,3);M=dct2(RGB1(:,:,3);M=(N-M)/Q;B=idct2(M(236:365,350:401);Y=mat2gray(B);imshow(Y);title('提取的水印圖像')運(yùn)行結(jié)果： 水印全過程（經(jīng)剪切檢測水?。?水印加入程序Q=input('請輸入放縮因子的值(建議小于1)：Q')subplot(3,3,1)RGB=imread

20、('MM','jpg');imshow(RGB);title('原始圖像');subplot(3,3,2)imshow(RGB(:,:,3);title('B分量');subplot(3,3,3)N=dct2(RGB(:,:,3);imshow(log(abs(N),),colorbar;title('B分量能量分布');subplot(3,3,4)I=imread('lena1','bmp');imshow(I);title('灰度水印圖像');subplot(3,

21、3,5)M=dct2(I);imshow(log(abs(M),),colorbar;title('水印能量分布');subplot(3,3,7)J=M(1:128,1:128);J(128:464,128:364)=0;J=rot90(J);J=rot90(J);J(465:800,365:600)=0;J=rot90(J);J=rot90(J);N=N+Q*J;K=idct2(N);RGB(:,:,3)=K;imshow(RGB);title('加入水印后圖像');subplot(3,3,8)I=imcrop(RGB,1 1 598 798);imshow(

22、I);subplot(3,3,9)%水印提取程序subplot(3,3,6)RGB1=imread('MM','jpg');J=RGB1(:,:,3);X=J(1:799,1:599);N=dct2(I(:,:,3);M=dct2(X);M=(N-M)/Q;B=idct2(M(337:464,237:364);Y=mat2gray(B);imshow(Y);title('經(jīng)放縮后提取的水印圖像')運(yùn)行結(jié)果： 水印全過程（經(jīng)空域壓縮檢測水?。┏绦蛟创a%水印加入程序Q=input('請輸入放縮因子的值(建議小于1)：Q')P=inp

23、ut('請輸入您所希望的圖像放縮系數(shù)值(建議取值不要小于0.5):P=')subplot(3,3,1)RGB=imread('浙江臺州學(xué)院','jpg');imshow(RGB);title('原始圖像');subplot(3,3,2)imshow(RGB(:,:,3);title('B分量');subplot(3,3,3)N=dct2(RGB(:,:,3);imshow(log(abs(N),),colorbar;title('B分量能量分布');subplot(3,3,4)I=imread(&#

24、39;lena1','bmp');imshow(I);title('灰度水印圖像');subplot(3,3,5)M=dct2(I);imshow(log(abs(M),),colorbar;title('水印能量分布');subplot(3,3,7)J=M(1:128,1:128);J(128:364,128:400)=0;J=rot90(J);J=rot90(J);J(365:600,401:750)=0;J=rot90(J);J=rot90(J);N=N+Q*J;K=idct2(N);RGB(:,:,3)=K;imshow(RGB)

25、;title('加入水印后圖像');subplot(3,3,8)I=imresize(RGB,P,'nearest');imshow(I);title('壓縮P倍圖像');subplot(3,3,9)J=imresize(I,1/P,'nearest');imshow(J);title('再放大P倍還原圖像')%水印提取程序subplot(3,3,6)RGB1=imread('浙江臺州學(xué)院','jpg');N=dct2(J(:,:,3);M=dct2(RGB1(:,:,3);M=(N

26、-M)/Q;B=idct2(M(236:365,350:401);Y=mat2gray(B);imshow(Y);title('經(jīng)放縮后提取的水印圖像')：運(yùn)行結(jié)果：4.2 圖像水印的dwt算法%以下是水印提取算法clear all;clc; %保存時間start_time=cputime;figure(1); %讀出原始圖像subplot(1,2,1);input=imread('2.jpg');imshow(input);title('原始圖像'); %讀出水印圖像subplot(1,2,2);watermarked_image=imread

27、('watermarked.bmp');imshow(watermarked_image,);title('水印圖像'); %三色分離input=double(input);watermarked_image=double(watermarked_image);inputr=input(:,:,1);watermarked_imager=watermarked_image(:,:,1);inputg=input(:,:,2);watermarked_imageg=watermarked_image(:,:,2);inputb=input(:,:,3);water

28、marked_imageb=watermarked_image(:,:,3); %水印圖像R的分解Cwr,Swr=WAVEDEC2(watermarked_imager,2,'haar'); %圖像R的分解Cr,Sr=WAVEDEC2(inputr,2,'haar'); %水印圖像G的分解Cwg,Swg=WAVEDEC2(watermarked_imageg,2,'haar'); %圖像R的分解Cg,Sg=WAVEDEC2(inputg,2,'haar'); %水印圖像B的分解Cwb,Swb=WAVEDEC2(watermarke

29、d_imageb,2,'haar'); %圖像B的分解Cb,Sb=WAVEDEC2(inputb,2,'haar'); %提取水印小波系數(shù)%提取水印R的小波系數(shù)r=0.06;for k=0:3 whr(k+1,:)=Cwr(1+size(Cwr,2)/4+k*size(Cwr,2)/16:. size(Cwr,2)/4+(k+1)*size(Cwr,2)/16)-. Cr(1+size(Cr,2)/4+k*size(Cr,2)/16:. size(Cr,2)/4+(k+1)*size(Cr,2)/16); wvr(k+1,:)=Cwr(1+size(Cwr,2)

30、/2+k*size(Cwr,2)/16:. size(Cwr,2)/2+(k+1)*size(Cwr,2)/16)-. Cr(1+size(Cr,2)/2+k*size(Cr,2)/16:. size(Cr,2)/2+(k+1)*size(Cr,2)/16); wdr(k+1,:)=Cwr(1+3*size(Cwr,2)/4+k*size(Cwr,2)/16:. 3*size(Cwr,2)/4+(k+1)*size(Cwr,2)/16)-. Cr(1+3*size(Cr,2)/4+k*size(Cr,2)/16:. 3*size(Cr,2)/4+(k+1)*size(Cr,2)/16);end

31、whr=(whr(1,:)+whr(2,:)+whr(3,:)+whr(4,:)/(4*r);wvr=(wvr(1,:)+wvr(2,:)+wvr(3,:)+wvr(4,:)/(4*r);wdr=(wdr(1,:)+wdr(2,:)+wdr(3,:)+wdr(4,:)/(4*r);war=(Cwr(1:size(Cwr,2)/16)-Cr(1:size(Cr,2)/16)/r; %提取水印小波系數(shù)%提取水印G的小波系數(shù)g=0.03;for k=0:3 whg(k+1,:)=Cwg(1+size(Cwg,2)/4+k*size(Cwg,2)/16:. size(Cwg,2)/4+(k+1)*si

32、ze(Cwg,2)/16)-. Cg(1+size(Cg,2)/4+k*size(Cg,2)/16:. size(Cg,2)/4+(k+1)*size(Cg,2)/16); wvg(k+1,:)=Cwg(1+size(Cwg,2)/2+k*size(Cwg,2)/16:. size(Cwg,2)/2+(k+1)*size(Cwg,2)/16)-. Cg(1+size(Cg,2)/2+k*size(Cg,2)/16:. size(Cg,2)/2+(k+1)*size(Cg,2)/16); wdg(k+1,:)=Cwg(1+3*size(Cwg,2)/4+k*size(Cwg,2)/16:. 3*

33、size(Cwg,2)/4+(k+1)*size(Cwg,2)/16)-. Cg(1+3*size(Cg,2)/4+k*size(Cg,2)/16:. 3*size(Cg,2)/4+(k+1)*size(Cg,2)/16);endwhg=(whg(1,:)+whg(2,:)+whg(3,:)+whg(4,:)/(4*g);wvg=(wvg(1,:)+wvg(2,:)+wvg(3,:)+wvg(4,:)/(4*g);wdg=(wdg(1,:)+wdg(2,:)+wdg(3,:)+wdg(4)/(4*g);wag=(Cwg(1:size(Cwg,2)/16)-Cg(1:size(Cg,2)/16)

34、/g; %提取水印小波系數(shù)%提取水印B的小波系數(shù)b=0.12;for k=0:3 whb(k+1,:)=Cwb(1+size(Cwb,2)/4+k*size(Cwb,2)/16:. size(Cwb,2)/4+(k+1)*size(Cwb,2)/16)-. Cb(1+size(Cb,2)/4+k*size(Cb,2)/16:. size(Cb,2)/4+(k+1)*size(Cb,2)/16); wvb(k+1,:)=Cwb(1+size(Cwb,2)/2+k*size(Cwb,2)/16:. size(Cwb,2)/2+(k+1)*size(Cwb,2)/16)-. Cb(1+size(Cb,2)/2+k*size(Cb,2)/16:. size(Cb,2)/2+(k+1)*size(Cb,2)/16); wdb(k+1,:)=Cwb(1+3*size(Cwb,2)/4+k*size(Cwb,2)/16:. 3*size(Cwb,2)/4+(k+1)*size(Cwb,2)/16)-. Cb(1+3*size(Cb,2)/4+k*size(Cb,2)/16:. 3*siz