%DEBLOCK   Remove blocking artifact from JPEG compressed images
%
% SYNOPSIS:
%  out = deblock(in)
%
% EXAMPLE:
%  a = readim
%  imwrite(uint8(a),'temp.jpg','jpg','Quality',10);
%  b = dip_image(imread('temp.jpg','jpg'))
%  c = deblock(b)
%  delete('temp.jpg')
%
% LITERATURE:
%  T.Q. Pham and L.J. van Vliet,  Blocking artifacts removal by
%   a hybrid filter method, submitted to ECCV'04
 
% (C) Copyright 1999-2003               Pattern Recognition Group
%     All rights reserved               Faculty of Applied Physics
%                                       Delft University of Technology
%                                       Lorentzweg 1
%                                       2628 CJ Delft
%                                       The Netherlands
%
% Tuan Pham, December 2003.

function out = deblock(in)

if nargin == 1
   if ischar(in) & strcmp(in,'DIP_GetParamList')
      out = struct('menu','none');
      return
   end
end


% jpeg_quality_score range=[0 10] with 10 being best quality. 
jqs = min(max(0,jpeg_quality_score(in)),10) % the score could get slightly>10

% the certainty image: 1 everywhere, <1 at block boundaries
block = newim([8 8])+jqs/10; block(1:6,1:6) = 1;
mask = repmat(block,ceil(size(in)/8));
c = gaussf(mask(0:size(in,1)-1,0:size(in,2)-1)); c = c/max(c);

% if input is perfect, s1d should be close to 0 to avoid blurring
s1d = (10-jqs)/7.5   % 7.5 is better than 10 for low bit-rate (connect lines better)
[p,q] = structuretensor(in,1,s1d,{'orientation','anisotropy'});
curv = curvature(in,'line',1,s1d);   % may skip this for speed
d = arcf(in,newimar(p,s1d*q,curv));

% for low freq region, use simple Gaussian smoothing, 
s2d = (10-jqs)/2     % sigma_2D range=[0 5]
e = gaussf(in,s2d);  

% normalized convolution to reduce blocking effect on gradmag
g = gaussf(sqrt(dx(in)^2 + dy(in)^2),s2d)/gaussf(c,s2d);    

% Gaussian error norm, Haglund SCIA'03 paper use sigmoid
w = exp(-g/median(g)/max(1,s2d-1));      
out = d*(1-w) + e*w ;

% compensate for interpolation blur caused by arcf
out = out - laplace(out)/2;




function score = jpeg_quality_score(img)

%========================================================================
%
%Copyright (c) 2002 The University of Texas at Austin
%All Rights Reserved.
% 
%This program is free software; you can redistribute it and/or modify
%it under the terms of the GNU General Public License as published by
%the Free Software Foundation; either version 2 of the License, or
%(at your option) any later version.
% 
%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.
% 
%The GNU Public License is available in the file LICENSE, or you
%can write to the Free Software Foundation, Inc., 59 Temple Place -
%Suite 330, Boston, MA 02111-1307, USA, or you can find it on the
%World Wide Web at http://www.fsf.org.
%
%Author  : Zhou Wang 
%Version : 1.0
% 
%The authors are with the Laboratory for Image and Video Engineering
%(LIVE), Department of Electrical and Computer Engineering, The
%University of Texas at Austin, Austin, TX.
%
%Kindly report any suggestions or corrections to zhouwang@ieee.org
%
%========================================================================
%
%This is an implementation of the algorithm for calculating the quality
%score of JPEG compressed images proposed by Zhou Wang, Hamid R. Sheikh
%and Alan C. Bovik. Please refer to the paper: Zhou Wang, Hamid R. Sheikh
%and Alan C. Bovik, "No-Reference Perceptual Quality Assessment of JPEG
%Compressed Images," submitted to IEEE International Conference on Image
%Processing, Sept. 2002.
%
%You can change this program as you like and use it anywhere, but please
%refer to its original source (cite our paper and our web page at
%http://anchovy.ece.utexas.edu/~zwang/research/nr_jpeg_quality/index.html).
%
%Input : A test 8bits/pixel grayscale image loaded in a 2-D array
%Output: A quality score of the image. The score typically has a value
%        between 1 and 10 (10 represents the best quality, 1 the worst).
%
%Usage:
%
%1. Load the image, for example
%
%   image = imread('testimage.jpg'); 
%
%2. Call this function to calculate the quality score:
%
%   Quality_Score = jpeg_quality_score(image)
%
%========================================================================

if (nargin > 1)
    score = -1;
    return;
end

x = double(img);

[M N] = size(x);
if (M < 16 | N < 16)
    score = -2;
    return;
end   

% Feature Extraction:

% 1. horizontal features

d_h = x(:, 2:N) - x(:, 1:(N-1));

B_h = mean2(abs(d_h(:, 8:8:8*(floor(N/8)-1))));

A_h = (8*mean2(abs(d_h)) - B_h)/7;

sig_h = sign(d_h);
left_sig = sig_h(:, 1:(N-2));
right_sig = sig_h(:, 2:(N-1));
Z_h = mean2((left_sig.*right_sig)<0);

% 2. vertical features

d_v = x(2:M, :) - x(1:(M-1), :);

B_v = mean2(abs(d_v(8:8:8*(floor(M/8)-1), :)));

A_v = (8*mean2(abs(d_v)) - B_v)/7;

sig_v = sign(d_v);
up_sig = sig_v(1:(M-2), :);
down_sig = sig_v(2:(M-1), :);
Z_v = mean2((up_sig.*down_sig)<0);

% 3. combined features

B = (B_h + B_v)/2;
A = (A_h + A_v)/2;
Z = (Z_h + Z_v)/2;

% Quality Prediction

alpha = -245.8909;
beta = 261.9373;
gamma1 = -239.8886;
gamma2 = 160.1664;
gamma3 = 64.2859;

score = alpha + beta*(B.^(gamma1/10000))*(A.^(gamma2/10000))*(Z.^(gamma3/10000));
