dct_and_idct_for_image.txt

dct and idct for image.matlab实现信息序列的嵌入和提取

%start by :  2007 09 24 14:30
%only show us how much can a dct and idct process change a image
%and what disturibute do the dct coefficient obey
%by wubaishan   version 1.1

%original_image:取成block的整数倍的原始图片
%one_block_dct：临时DCT变换所取的块的DCT系数
%DCT2:按2维存储的（M×N）×64 的系数矩阵
%DCT3：按3维存储的M*N×64的zig－zag排列过的系数矩阵 size:row_BlockNum*column_BlockNum
%DCT4:按4维存储的M×N × 8×8 的系数矩阵 size:row_BlockNum*column_BlockNum*8*8
%DCT_:处理过的DCT     size:row_BlockNum*column_BlockNum
%
%
%
%
%




clc
clear all
close all

%parameters
quantity_scale=0.1;
T_psnr=40;
K=16;
test=0.97;
%end parameters

%constant define
BlockWidth=8;
quantity_width=4;
bm=255;
%end constant define

%constant matrix
permute_matrix_form1=[
   0  1  5  6  14 15 27 28;
   2  4  7  13 16 26 29 42;
   3  8  12 17 25 30 41 43;
   9  11 18 24 31 40 44 53;
   10 19 23 32 39 45 52 54;
   20 22 33 38 46 51 55 60;
   21 34 37 47 50 56 59 61;
   35 36 48 49 57 58 62 63;
   ];
                    
quantity_standrad_matrix=[
    16  11  10  16  24  40  51  61 ;
    12  12  14  19  26  58  60  55 ;
    14  13  16  24  40  57  69  56 ; 
    14  17  22  29  51  87  80  62 ;
    18  22  37  56  68  109 103 77 ;
    24  35  55  64  81  104 113 92 ;
    49  64  78  87  103 121 120 101;
    72  92  95  98  112 100 103 99 ;
];
%end constant matrix

%calculate constant parameters
quantity_scale_matrix=quantity_scale*quantity_standrad_matrix;%used to decise the quantity steps

permute_matrix_form2=matrix_change(permute_matrix_form1);%used to achieve zig-zag order(from 2D to 1D)
%let quantity matrix into a line
for quantity_i=1:8
    for quantity_j=1:8
        quantity_scale_line((quantity_i-1)*8+quantity_j)=quantity_scale_matrix(quantity_i,quantity_j);
    end
end

%get gray image
original_image1=imread('city7.bmp');     %read image matrix
[high,width,is_color]=size(original_image1);%m is high ,n is width
if is_color==3                              %if it's colorful image,change it to gray form
    original_image1=rgb2gray(original_image1);
end

%get image as a multiply of 8*8
column_BlockNum=fix(width/BlockWidth)   %how many blocks can this original image divide to
row_BlockNum=fix(high/BlockWidth)       %how many blocks can this original image divide to
original_image=original_image1(1:row_BlockNum*BlockWidth,1:column_BlockNum*BlockWidth); %get a image of interal multiply of BlockWidth
original_image=double(original_image);  %matrix get from imread is the size of uint8

%calculate constant parameters
    
mark_strength=test*8*bm*10^(-T_psnr/20)/(K)^0.5      %insert strength

%get gaussian distribution of 0 and 1
watermark_1=randn(1,row_BlockNum*column_BlockNum);
for watermark_inf_i=1:row_BlockNum*column_BlockNum
    if watermark_1(watermark_inf_i)>0
        watermark_inf(watermark_inf_i)=1;
    else
        watermark_inf(watermark_inf_i)=0;
    end
end 
   
%end calculate constant parameters


%dct change for each block
for i=0:row_BlockNum-1        %row point            %outer block loop
    for j=0:column_BlockNum-1 %column point1         %outer block loop
        
        %get pixel matrix data (BlockWidth by BlockWidth)
        for k=1:BlockWidth    %column point2
            pixel_matrix(k,:)=original_image(i*BlockWidth+k,j*BlockWidth+1:(j+1)*BlockWidth);
        end   
        %end get pixel matrix data (BlockWidth by BlockWidth)
        
        %dct change
        one_block_dct=dct2(pixel_matrix);
        %end dct change
              
        %lie those coefficients into zig-zag line 
        for l=1:BlockWidth^2
            zig_zag_line_dct(l)=one_block_dct(permute_matrix_form2(l));
        end
        %end lie those coefficients into zig-zag line
        
        %storage each block dct coefficient to the big DCT matrix(3D),DCT1(2D)
        DCT2(i*row_BlockNum+j+1,:)=zig_zag_line_dct(:); %DCT1 is only for display and it's 2D
        DCT3(i+1,j+1,:)=zig_zag_line_dct(:);
        DCT4(i+1,j+1,:,:)=one_block_dct(:,:);
        %end storage each block dct coefficient to the big DCT matrix(3D),DCT1(2D)
        
    end     %end :for j=0:column_BLockNum-1 %column point         %outer block loop
end    %for i=0:row_BlockNum-1        %row point          







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%this time,we mark the inf directly in the 10 bigger coefficient(expect the DC component)
%mark_strength=0;

%initialization
DCT_=DCT3;

for mark_i=1:row_BlockNum
    for mark_j=1:column_BlockNum
        %if watermark_inf((mark_i-1)*column_BlockNum+mark_j)==1
            DCT_(mark_i,mark_j,1:1+K)=DCT3(mark_i,mark_j,1:1+K)+mark_strength;
        %end
    end
end
%end this time,we mark the inf directly in the DC coefficient
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DCT_=DCT3;%if no noise was additive


%idct process
for idct_i=1:row_BlockNum               %outer block loop (Row loop)
    for idct_j=1:column_BlockNum        %outer block loop (colume loop)
        
        %inv zig-zag transform,coefficients was get form DCT_ matrix
        for permute_i=1:BlockWidth
            for permute_j=1:BlockWidth
                permute_coefficient(permute_i,permute_j)=DCT_(idct_i,idct_j,permute_matrix_form1((permute_i-1)*BlockWidth+permute_j)+1);
            end
        end
        %end inv zig-zag transform
        
        %block idct change
        watermark_block_image=idct2(permute_coefficient);
        %end block idct change
                
        %storage the block_image to watermark_image matrix
        for together_i=1:BlockWidth
            watermark_image((idct_i-1)*BlockWidth+together_i,(idct_j-1)*BlockWidth+1:(idct_j-1)*BlockWidth+BlockWidth)=watermark_block_image(together_i,:);
        end
        %end storage the block_image to watermark_image matrix
        
    end
end


%simulation results
%original_image=fix(original_image);
%watermark_image=fix(watermark_image);
image_dif=original_image-watermark_image;
%error_num=nnz(fix(image_dif))
%error_rate=error_num/(row_BlockNum*BlockWidth*column_BlockNum*BlockWidth)


image_dif_mse=mse(double(image_dif));
original_image_mse=mse(double(original_image));
SNR=20*log10(original_image_mse/image_dif_mse)
PSNR=20*log10(bm)-10*log10(image_dif_mse)









%figure(1)
%original_image_show=uint8(original_image);
%imshow(original_image_show)
%figure(2)
%watermark_image_show=uint8(watermark_image);
%imshow(watermark_image_show)