first.m

dct based visible watermarking algorithm

clc;
clear all;
close all;
I=imread('original1.jpg'); %Get the input image
W=imread('wat.jpg');%watermark image

I=rgb2gray(I);      %Convert to grayscale image
W=rgb2gray(W);

I=im2double(I);
W=im2double(W);

figure,imshow(I),title('Original Image');
[p q] = size(I);

figure,imshow(W),title('Watermark Image');
[r s] = size(W);
n1=p*q/64; %no of blocks in original image
n2=r*s/64; %no of blocks in watermark image
l=p/8;
m=q/8;

%Divide the images I and W into 8*8 blocks

c=0;
    for b=1:m
     for a=1:l
      for j=1:8
       for k=1:8
        Block(a+l*c).X1(j,k)=I(j+8*(a-1),k+8*(b-1));
       end
      end
     end
     c = c+1;
    end

    
l=r/8;
m=s/8;
c=0;
 for b=1:m
  for a=1:l
   for j=1:8
    for k=1:8
     Block(a+l*c).Y1(j,k)=W(j+8*(a-1),k+8*(b-1));
    end
   end
  end
 c = c+1;
 end


 % Calculate DCT coefficients for both images

 for i=1:n1 
     Block(i).X2=dct2(Block(i).X1);
 end


 for i=1:n2
     Block(i).Y2=dct2(Block(i).Y1);
 end

 
 % Calculate Normalized mean grey value (nmgv) for each Block of image I

 MaxVal = Block(1).X1(1,1);  %Assign Maximum value to first block,first pixel that is the DC DCT coefficient

for i=2:n1
  if(Block(i).X1(1,1)>MaxVal)
      MaxVal = Block(i).X1(1,1);
  end      
end

%Normalized Mean Gray value os calculated for each block (using eq 4 of the paper)
for i = 1:n1
    Block(i).Norm = Block(i).X1(1,1) / MaxVal;
end

Sum = 0;

for i = 1:n1
   Sum = Sum + Block(i).X1(1,1);
end

nmgv = Sum / n1; %nmgv stores the normalized image mean gray value (eq 5)

%Calculate Normalized variance for each block of I
Sum = 0;

for i = 1:n1
    Sum = Sum + Block(i).X2;
end

MeanDct = Sum/n1;

for i = 1:n1
    Block(i).VarDct = ((Block(i).X2 - MeanDct)* (Block(i).X2 - MeanDct))/64 ;
    Block(i).LogVarDct = log(Block(i).VarDct) ;
end

MaxVar = Block(1).LogVarDct;

for i = 1:n1
     if (Block(i).LogVarDct > MaxVar)
          MaxVar = Block(i).LogVarDct ;
     end
end

for i = 1:n1
    Block(i).NorVar = Block(i).LogVarDct/MaxVar;
end

%Identification of Edge Blocks using Sobel method

BW = edge(I,'sobel');
figure,imshow(BW);
       
c = 1;

% Calculating Scale Factor and Embedding Factor of each block
for i =1:n1
     Block(i).ScaleFact = Block(i).LogVarDct * exp ( - ((Block(i).Norm - nmgv)*(Block(i).Norm - nmgv)));
end

for i =1:n1
    a=exp ( - ((Block(i).Norm - nmgv)*(Block(i).Norm - nmgv)));
    y=1-a;
    z=Block(i).LogVarDct.^(-1);
     Block(i).EmbedFact = y*z;
end

%Find maximum value of scaling factor
MaxScaleFact = Block(1).ScaleFact;

for i = 1:n1
    if (MaxScaleFact < Block(i).ScaleFact)
     MaxScaleFact = Block(i).ScaleFact;
    end
end

%Find Minimum value of embedding factor
MinEmbedFact = Block(1).EmbedFact;

for i = 1:n1
    if(MinEmbedFact > Block(i).EmbedFact)
        MinEmbedFact = Block(i).EmbedFact;
    end
end

%Divide the egde image (BW) into 8*8 blocks ( to identify the edge blocks)
l=p/8;
m=q/8;
c=0;

for b=1:m
   for a=1:l
      for j=1:8
       for k=1:8
        Block(a+l*c).X5(j,k)=BW(j+8*(a-1),k+8*(b-1));
       end
      end
     end
     c = c+1;
    end
  

  % Assign previously calculated Scaling factors and Embedding Factors for each block
    for i = 1:n1
        for j = 1:8
            for k = 1:8
                if (Block(i).X5(j,k) == 1) %since the image is binary, the gray value will be one for the edge block and zero otherwise
                    Block(i).ScaleFact = MaxScaleFact ;
                    Block(i).EmbedFact = MinEmbedFact;
                  
                else
                    Block(i).ScaleFact = Block(i).LogVarDct * exp ( - ((Block(i).Norm - nmgv)*(Block(i).Norm - nmgv)));
                    a=exp ( - ((Block(i).Norm - nmgv)*(Block(i).Norm - nmgv)));
                    y=1-a;
                    z=Block(i).LogVarDct^(-1);
                    Block(i).EmbedFact = y*z;
                   
                end
            end
        end
    end

%Calculate DCT coefficients of the watermarked image

for i=1:n1
    Block(i).NewDct = Block(i).ScaleFact*Block(i).X2 + Block(i).EmbedFact*Block(i).Y2;
%    disp(Block(i).NewDct);
end

%Calculate Inverse DCT of watermarked image
for i=1:n1
    Block(i).X3=idct2(Block(i).NewDct);
end  

%Unite the image (Reverse of division)

c=0;
l=p/8;
m=q/8;
for b=1:m
    for a=1:l
      for j=1:8
       for k=1:8
        Wimg(j+8*(a-1),k+8*(b-1))=Block(a+l*c).X3(j,k);
       end
      end
     end
     c = c+1;
end
figure,imshow(Wimg),title('watermarked Image');