1. 程式人生 > >基於DCT變換(變換域)實現資訊隱藏(數字水印)原理詳解及matlab實現

基於DCT變換(變換域)實現資訊隱藏(數字水印)原理詳解及matlab實現

主要就是實現了數字水印的嵌入提取和在不同攻擊如旋轉剪下噪聲等下的提取效果差異的比對

1 DCT變換的原理

2  DCT變換的特點

    在基於DCT的變換編碼中,影象是先經分塊(8×8或16×16)後再經DCT,這種變換是區域性的,只反映了影象某一部分的資訊。當然也可以對整幅影象的特點,但是運算速度比分塊DCT要慢。影象經DCT後,得到的DCT影象有三個特點:  

    (1). 係數值全部集中到0值附近(從直方圖統計的意義上),動態範圍很小,這說明用較小的量化位元數即可表示DCT係數; 

    (2). DCT變換後圖像能量集中在影象的低頻部分,即DCT影象中不為零的係數大部分集中在一起(左上角),因此編碼效率很高。 

    (3). 沒有保留原影象塊的精細結構,從中反映不了原影象塊的邊緣、輪廓等資訊,這一特點是由DCT缺乏時局域性造成的。

3 宿主影象的DCT變換

對於N×N大小的256灰度級的宿主影象I進行N×N二維離散餘弦變換(DCT)。以ZigZag方式對於DCT變換後的影象頻率係數重新排列成一維向量Y={y1, y2,…yN×N}.  並取出序列中第L+1到L+M的中頻係數部分,得到YL={ YL+1, YL+2,…, YL+M}。

4 DCT數字水印嵌入演算法流程


5.DCT數字水印提取演算法流程


matlab程式碼實現:

DCT主函式

%------------------------------------------------------------------%
%     基於DCT變換的資訊隱藏(數字水印)     %%                                                                %%                                                                %
%-----------------------------------------------------------=------%
clear ;
clc;
      
%-----------------------讀入影象-------------------------------------%
markbefore=imread('sss.bmp');
markbefore2=rgb2gray(markbefore);
mark=im2bw(markbefore2);    %使水印影象變為二值圖
figure(1);      %開啟視窗
subplot(2,3,1);    %該視窗內的影象可以有兩行三列
imshow(mark),title('水印影象');   %顯示水印影象
[rm,cm]=size(mark);   %計算水印影象的長寬

cover_image=imread('carrier_image.bmp');
cover_image=rgb2gray(cover_image);
subplot(2,3,2),imshow(cover_image,[]),title('載體影象'); %[]表示顯示時灰度範圍為image上的灰度最小值到最大值

before=blkproc(cover_image,[8 8],'dct2');   %將載體影象的灰度層分為8×8的小塊,每一塊內做二維DCT變換,結果記入矩陣before

I=mark;
alpha=30;     %尺度因子,控制水印新增的強度,決定了頻域係數被修改的幅度
k1=randn(1,8);  %產生兩個不同的隨機序列
k2=randn(1,8);
after=before;   %初始化載入水印的結果矩陣
for i=1:rm          %在中頻段嵌入水印
    for j=1:cm
        x=(i-1)*8;
        y=(j-1)*8;
        if mark(i,j)==1
            k=k1;
        else
            k=k2;
        end;
        after(x+1,y+8)=before(x+1,y+8)+alpha*k(1);
        after(x+2,y+7)=before(x+2,y+7)+alpha*k(2);
        after(x+3,y+6)=before(x+3,y+6)+alpha*k(3);
        after(x+4,y+5)=before(x+4,y+5)+alpha*k(4);
        after(x+5,y+4)=before(x+5,y+4)+alpha*k(5);
        after(x+6,y+3)=before(x+6,y+3)+alpha*k(6);
        after(x+7,y+2)=before(x+7,y+2)+alpha*k(7);
        after(x+8,y+1)=before(x+8,y+1)+alpha*k(8);
    end;
end;
result=blkproc(after,[8 8],'idct2');    %將經處理的影象分為8×8的小塊,每一塊內做二維DCT逆變換
result = uint8(result);
imwrite(result,'markresule.bmp','bmp');      %儲存新增水印後的影象
subplot(2,3,3),imshow(result,[]),title('嵌入水印的影象');    %顯示新增水印後的影象

%定義一個空空間來儲存提取的水印
disp('請選擇對影象的攻擊方式:');
disp('1.新增白噪聲');
disp('2.高斯低通濾波處理');
disp('3.對影象進行部分剪下');
disp('4.將影象旋轉十度');
disp('5.將影象壓縮處理');
disp('6.新增椒鹽噪聲');
disp('7.不處理影象,直接顯示提取水印');
disp('輸入其它數字則直接顯示提取水印');
choice=input('請輸入選擇:');
figure(1);
switch choice        %讀入輸入的選擇  withmark為等待提取水印的影象
    case 1
        result_1=result;
        withmark=imnoise(result_1,'gaussian',0.02); %加入椒鹽躁聲
        subplot(2,3,4);
        imshow(withmark,[]);
        title('加入高斯白噪聲後的影象');     %顯示加了椒鹽噪聲的影象
    case 2
         result_2=result;
         H=fspecial('gaussian',[4,4],0.2); 
         result_2=imfilter(result_2,H); 
         subplot(2,3,4); 
         imshow(result_2,[]); 
         title('高斯低通濾波後圖像'); 
         withmark=result_2;
    case 3
        result_3=result;
        result_3(1:64,1:400)=512;   %使影象上方被剪裁
        subplot(2,3,4);
        imshow(result_3);
        title('上方剪下後圖像');
        withmark=result_3;
    case 4
        result_4=imrotate(result,10,'bilinear','crop');   %最鄰近線性插值演算法旋轉10度
        subplot(2,3,4);
        imshow(result_4);
        title('旋轉10度後圖像');
        withmark=result_4;
    case 5
        result_5 = result; 
        result_5=im2double(result_5); 
        cnum=10; 
        dctm=dctmtx(8); 
        P1=dctm; 
        P2=dctm.'; 
        imageDCT=blkproc(result_5,[8,8],'P1*x*P2',dctm,dctm.'); 
        DCTvar=im2col(imageDCT,[8,8],'distinct').'; 
        n=size(DCTvar,1); 
        DCTvar=(sum(DCTvar.*DCTvar)-(sum(DCTvar)/n).^2)/n; 
        [dum,order]=sort(DCTvar); 
        cnum=64-cnum; 
        mask=ones(8,8); 
        mask(order(1:cnum))=zeros(1,cnum); 
        im88=zeros(9,9); 
        im88(1:8,1:8)=mask; 
        im128128=kron(im88(1:8,1:8),ones(16)); 
        dctm=dctmtx(8); 
        P1=dctm.'; 
        P2=mask(1:8,1:8); 
        P3=dctm; 
        result_5=blkproc(imageDCT,[8,8],'P1*(x.*P2)*P3',dctm.',mask(1:8,1:8),dctm); 
        WImage5cl=mat2gray(result_5); 
        subplot(2,3,4); 
        imshow(WImage5cl); 
        title('經JPEG壓縮後圖像'); 
        withmark=WImage5cl;
    case 6
        result_6=result;
        withmark=imnoise(result_6,'salt & pepper',0.02); %加入椒鹽躁聲
        subplot(2,3,4);
        imshow(withmark,[]);
        title('加入椒鹽噪聲後的影象');     %顯示加了椒鹽噪聲的影象
    case 7
        subplot(2,3,4);
        imshow(result,[]);
        title('未受攻擊的水印影象');
        withmark=result;
    otherwise
        disp('輸入數字選擇無效,影象未受攻擊,直接提取水印');
        subplot(2,3,4);
        imshow(result,[]);
        title('未受攻擊的水印影象');
        withmark=result;
end

%------------------------水印提取-----------------------------%
%
after_2=blkproc(withmark,[8,8],'dct2');   %此步開始提取水印,將灰度層分塊進行DCT變換
p=zeros(1,8);        %初始化提取數值用的矩陣
mark_2 = zeros(rm,cm);
for i=1:rm
    for j=1:cm
        x=(i-1)*8;y=(j-1)*8;
        p(1)=after_2(x+1,y+8);         %將之前改變過數值的點的數值提取出來
        p(2)=after_2(x+2,y+7);
        p(3)=after_2(x+3,y+6);
        p(4)=after_2(x+4,y+5);
        p(5)=after_2(x+5,y+4);
        p(6)=after_2(x+6,y+3);
        p(7)=after_2(x+7,y+2);
        p(8)=after_2(x+8,y+1);
        if corr2(p,k1)>corr2(p,k2)  %corr2計算兩個矩陣的相似度,越接近1相似度越大
            mark_2(i,j)=1;              %比較提取出來的數值與隨機頻率k1和k2的相似度,還原水印圖樣
        else
            mark_2(i,j)=0;
        end
    end
end
subplot(2,3,5);
mark_2 = uint8(mark_2);
imshow(mark_2,[]),title('提取出的水印');
subplot(2,3,6);
imshow(mark),title('原嵌入水印');
NC=correlation(mark_2,mark);  
disp('原水印影象與提取水印影象互相關係數:')
disp(NC);
correlation函式:
function N=correlation(mark_get,mark_prime) 
mark_get=double(mark_get); 
mark_prime=double(mark_prime); 
if size(mark_get)~=size(mark_prime) 
    error('Input vectors must  be the same size!') 
else 
    [m,n]=size(mark_get); 
    fenzi=0; 
    fenmu=0; 
    for i=1:m 
        for j=1:n 
            fenzi=fenzi+mark_get(i,j)*mark_prime(i,j); 
            fenmu=fenmu+mark_prime(i,j)*mark_prime(i,j); 
        end 
    end 
N=min(fenzi/fenmu,fenmu/fenzi); 
end
其中的一個效果圖