% layer_model.m
%	check the layer model for DCT2 coefficients of images
%
% Written: Hao Fang
% Created: Aug. 2012

close all;
clear;
clc;

path(path, '../Tools/ZIGZAG_SCAN');
path(path, '../../Images');

img = imread('boat.512.tiff'); %r0=32; alpha=0.2
% img = imread('cameraman.tif'); %r0=14; alpha=0.4
% img = imread('Lena.tiff'); %r0=30; alpha=0.15
% img = rgb2gray(img);
% img = imread('Peppers.tiff'); %r0=38; alpha=0.15
% img = rgb2gray(img);

alpha = 0.15;

[M N] = size(img);
img_dct2 = dct2(img);
img_dct2_zigzag = zigzag(img_dct2);

threshold = 10e2;
r0 = 50;
% r = r0:M;
r = 1:r0;
ind_r = (r + 1) .*r ./2;
p = zeros(1, length(r));

img_dct2_zigzag_t = abs(img_dct2_zigzag) >= threshold;
sparsity = sum(img_dct2_zigzag_t) / (M*N);

num_of_large = zeros(1, length(r));
for ii = 1:length(r)
    num_of_large(ii) = sum(img_dct2_zigzag_t(ind_r(ii)-r(ii)+1:ind_r(ii)));
    p(ii) = num_of_large(ii) / r(ii);
end

r1 = 3;
r2 = 32;
y = exp(-alpha*(r-1));
y(1:r1) = 1;
y(r2+1:end) = 0;
figure('Position',[50 50 800 500]);
hold all;
plot(p, '-s','LineWidth', 2);
hold all;
plot(y, '-o','LineWidth', 2);

set(gca, 'FontSize', 14, ...
    'FontWeight', 'Bold');

xlabel('Layer Index m', 'FontSize', 20, 'FontWeight', 'Bold');
ylabel('p_m', 'FontSize', 20, 'FontWeight', 'Bold');

h_legend = legend('real', 'model');
set(h_legend, 'FontSize', 15);
set(h_legend, 'FontWeight', 'Bold');

