function [hidden_image, no_bits_hidden]=zeroDivJPEGstego(filename,per,design_quality, ...
    output_quality,threshold,noused,dseedint,bseedint,pseedint)

%% ################################################
%% ################################################
%% Author: Kaushal M. Solanki
%% Dept. of ECE, UCSB
%% Code developed in Sept 2005
%% solanki@ece.ucsb.edu
%% ################################################
%% ################################################
%% example inputs--
%% filename='C:\Test Images\PEPPERS.TIF';
%% per=30; %% per denotes the percentage of the terms actually used for
%% hiding
%% design_quality=50; %% system is designed for an input quality factor 
%% output_quality=75; %% we can resist compression when JPEG compression
%% uses the output quality factor > input quality factor
%% threshold=30; %% we hide data in AC DCT terms of magnitude <= threshold
%% noused=19; %% The first 19 AC DCT coefficients are used per 8x8 block,
%% that occur  in zigzag scan order
%% dseedint=1000; %% seed for generating the input dither sequence
%% bseedint=2000; %% seed for generating the actual bitstream
%% pseedint=3000; %%initial seed for random permutation vector to decide
%% data hiding locations


% JPEG steganography
% Input image is a jpeg image at a particular quality factor.
% The goal is to output another JPEG image which resembles the original
% image statistically (at that quality factor), and would still survive a design quality factor
% recompression attack, provided that the design quality factor is smaller than output
% quality factor

% Define JPEG quantization matrix
designJQM = getjqm(design_quality);
outputJQM = getjqm(output_quality);
dmin=1; % parameter indicating the binwidth used in constructing the PMF from continuous DCT distribution

% get the input image 
lena = imread(filename);

%%convert RGB to YUV.. embed data only in Y component
if size(lena,3)==3
    yuv = rgb2ntsc(lena);
    lena = uint8(255*yuv(:,:,1));
end

rows = size(lena,1);
cols = size(lena,2);

% number of 8x8 blocks in the image
noblocks = floor(rows/8)*floor(cols/8);

% Available coefficients for hiding and compensation
avcoef = noblocks*noused;

% payload length: number of bits used only in hiding
nobits = floor((per*avcoef)/100);

% generate the used matrix, which is used to select the AC DCT terms being
% used for hiding
usedt = genused(noused);
used = usedt(1:8,1:8);
clear usedt

% Create the dither sequence
dseed = dseedint; %%3142; % to be shared by encoder and decoder
rand('state',dseed);
ditherlevel = 1*dmin;
dseq = ditherlevel*(rand(nobits,1)-dmin/2);

% set a new seed for random bitstream: bit denotes the input binary message
% sequence
bseed = bseedint; %52105;
rand('state',bseed)
% bit = randint(nobits,1); % generate the bitstream
bit = round(rand(nobits,1)); 

% generate the permutation for locations of hiding vs. those for
% compensation 
vsel = zeros(avcoef,1);
vsel(1:nobits)=1;

% seed for the random permutation. This is for selecting which coefficients
% to hide and which to use for compensation.
pseed = pseedint; %%5190;
rand('state',pseed)
prm = randperm(avcoef);
sel=zeros(avcoef,1);
for k=1:avcoef
    sel(k)=vsel(prm(k));
end
clear vsel

% embedding process: hlena is the hidden image
hlena = lena;
hdctimg = zeros(rows,cols);
count = 0; % counter for locations for hiding and compensation
hcount =0; % count for hiding coefficients
ccount =0; % count for compensation coefficients

hStream = zeros(nobits,1);      % hiding stream
cStream = zeros(avcoef-nobits,1);    % compensation stream
stream = zeros(avcoef,1); % avilable stream
emb = zeros(avcoef,1);
hdistortion=0; % hiding induced distortion
abc=0;

% An important parameter: the threshold 
%% hiding occurs in DCT coefficients with magnitude <= threshold
thres = threshold; %%dmin*30;


%% most distributions encountered in practice, such as Gaussian or
%% Laplacian, have low prob tails and it is possible to avoid embedding in
%% low prob. region by using a threshold , i.e. the encoder would not embed
%% in the host symbols with absolute values greater than a pre determined 
%% threshold.. the decoder shares this threshold value which then uses the
%% same criteria for decoding..

for i=1:floor(rows/8)
    for j=1:floor(cols/8)
        block = lena((i-1)*8+1:i*8,(j-1)*8+1:j*8); 
        dctblock = dct2(block);
        %jplain=(dctblock./outputJQM);
        jblock = round(dctblock./outputJQM);
        hjblock = jblock;
        for m=1:8
            for n=1:8
                if used(m,n) == 1
                    p = jblock(m,n);
                    count=count+1;                   
                    stream(count)=p; 
                                       
                    if sel(count)==1       
                        hcount = hcount+1;
                        sc = designJQM(m,n)/outputJQM(m,n);
                        ds = sc*dseq(hcount);  
                        sdmin = sc*dmin;
                        
                        %% Those coeff whose magnitude greater than an
                        %%integer threshold are not used for hiding..
                        %%as they are in low prob. regions
                        if (p<=round(thres+ds))&&(p>=round(-thres+ds)) 
                            abc=abc+1;
                            
                            %%convert to odd no to hide 1 and to even no to
                            %%hide 0..we cannot use  SEC as p is whole no
                            %%and not fraction
                            if bit(hcount) == 1
                                q = p+sdmin-mod(p-ds,2*sdmin);
                            else
                                q = p+sdmin-mod(p+sdmin-ds,2*sdmin);
                            end
                            % Trying
                            q = round(q);
                            
                            % End trying
                            
                            hjblock(m,n) = q;
                            emb(count)=q;
                            hStream(hcount)=q;
                            %hbincount(q)=hbincount(q)+1;
                            hdistortion = hdistortion + (hjblock(m,n)-jblock(m,n))^2;
                           
                        else
                            
                            %% for low prob regions, we do not hide data
                            %% the idea is to match original histogram
                            %% precisely in high prob region
                            %% since low prob region is hard to
                            %% compensate,we just avoid embedding in that
                            %% region.. this way, the low prob region need
                            %% not be restored and we can get good
                            %% embedding rates by matching histograms in
                            %% high prob regions
                            hStream(hcount)=p;
                            emb(count)=p;
                        end                       
                    else   
                        emb(count)=p;
                        ccount = ccount+1;
                        cStream(ccount)=p;
                    end
                end
            end
        end
        
                        
        % store the current version of hjblock
        hdctimg((i-1)*8+1:i*8,(j-1)*8+1:j*8)=hjblock;
        hdctblock=hjblock.*outputJQM;
        hblock = uint8(round(idct2(hdctblock)));
        hlena((i-1)*8+1:i*8,(j-1)*8+1:j*8) = hblock;
    end
end
     
no_bits_hidden = abc;
disp('number of bits hidden..');



minStr = min(min(stream),min(emb));
maxStr = max(max(stream),max(emb));

x=minStr:maxStr;

x_display = -thres:thres; %%-30:30;

Hist2Comp = hist(stream,x)-hist(hStream,x);

NewStream = cStream; % The new variable for modified cStream

[scomp cindex] = sort(cStream);

maxdiff = thres; %%30;

mxd=0;
nco=0;

iComp=1;
last = length(x);
% Compensation procedure
binVal=minStr-1;% value in the current bin
maxComp = avcoef-nobits;% number of compensation coeffs

%iComp_broken = 0;

for i = 1:last
    binVal = binVal+1;
    for m = 1:Hist2Comp(i)
        if iComp > maxComp
            %iComp_broken = 1;
            break;
        end
        if (scomp(iComp) - binVal)^2 > 0.0002 % compensate
            if abs(binVal-scomp(iComp))<=maxdiff
                NewStream(cindex(iComp))=binVal;
            else
                mxd = mxd+1;
            end
            iComp = iComp+1; 
        else
            nco = nco+1;
            iComp = iComp+1; % go to the next compensation coefficient
        end
    end
    if iComp > maxComp
        %iComp_broken = 1;
        break;
    end   
end

ccount=0;
count=0;
for i=1:floor(rows/8)
    for j=1:floor(cols/8)
        hjblock = hdctimg((i-1)*8+1:i*8,(j-1)*8+1:j*8);
        for m=1:8
            for n=1:8
                if used(m,n) == 1
                    count=count+1;
                    if sel(count)~=1       
                        ccount=ccount+1;
                        hjblock(m,n)=NewStream(ccount);
                    end
                end
            end
        end
        hdctblock=hjblock.*outputJQM;
        % hdctblock=hjblock;
        hblock = uint8(round(idct2(hdctblock)));
        hlena((i-1)*8+1:i*8,(j-1)*8+1:j*8) = hblock;
    end
end

hidden_image = uint8(hlena);
figure,imshow(lena); title('Original image before steganography');
figure,imshow(hidden_image); title('Hidden image after steganography');

figure;hist(stream,x_display);title('Original histogram')
figure;hist(emb,x_display);title('Hidden but uncompensated histogram')
figure;bar(x_display,hist(stream,x_display)-hist(hStream,x_display));title('Target histogram for compensation Stream');
figure;bar(x_display,hist(NewStream,x_display)+hist(hStream,x_display)-hist(stream,x_display));title('difference')
figure;bar(x_display,hist(NewStream,x_display)+hist(hStream,x_display));title('final histogram')