deblock.m

jpeg的matlab代码

% deblock - smooth blocking artifacts
%
% SYNOPSIS
%	deblock(X)
%	deblock(X, ...)
%
% DESCRIPTION
% 	deblock performs a smoothing procedure in an effort to reduce
% 	blocking artifacts (the telltale discontinuities between blocks
%	which often follow aggressive quantizing). For a given block, a
%	polynomial is fitted using the DC coefficients from nearest 
%	neighbors. The k (default 5) lowest-frequency zero AC coefficients 
%	are replaced by values from the polynomial. Replaced values are 
%	clamped so that (in magnitude) they do not exceed values which 
%	would not have quantized to zero.  If the quantizer Q is not 
%	given, then the global QMAT is used.
%
%	The procedure is applied after dequantizing the quantized
%	transformed image (e.g., X = dequant(quant(dct(image), Q)) ).
%	The basic procedure is described in Pennebaker and Mitchell.
%
% BUGS
%	It's slow.

% Copyright (C) 1995-1997 Darrel Hankerson and Greg A. Harris
%
% This file is part of the JPEGtool collection of scripts for Octave
% and Matlab; see http://www.dms.auburn.edu/compression
%
% The JPEGtool collection 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, or
% (at your option) any later version.
%
% The collection 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.
%
% You should have received a copy of the GNU General Public License
% along with this package; see the file COPYING.  If not, write to the
% Free Software Foundation, 59 Temple Place - Suite 330, Boston, MA
% 02111-1307, USA.

function X = deblock(X, arg2, arg3)

global QMAT

if exist('arg2')
  if (length(arg2) > 1) Q = arg2; else k = arg2; end
  if exist('arg3')
    if (length(arg3) > 1) Q = arg3; else k = arg3; end
  end
end
    
if ~exist('Q')
  Q = QMAT;
end

% The number of entries to smooth can be specified (default is k=5).
if ~exist('k')
  k = 5;
end

N = size(Q, 1);

% Let 
%   p(v,u) = A1 u^2 v^2 + A2 u^2 v + A3 u v^2 + A4 u^2 + A5 u v
%            + A6 y^2 + A7 u + A8 v + A9
% be a surface over a 3x3 grid (where each element is an NxN block).
% Choose the coefficients A1--A9 so that the average value of p over each 
% element is the corresponding DC coefficient.
%
% Write in matrix form
%   MA = D
% and then find the matrix inverse of M (since we must solve for many D).
% Using an LU decomposition (rather than the inverse of M) may be preferable.

M = zeros(9);
for v = 1:9
  y = floor((v-1)/3)*N + 1; y = y:y+N-1;
  x = rem(v-1, 3)*N + 1; x = x:x+N-1;
  x1 = x * ones(N,1); y1 = y * ones(N,1);
  x2 = x * x'; y2 = y * y';
  M(v,:) = [x2*y2 x2*y1 x1*y2 N*x2 x1*y1 N*y2 N*x1 N*y1 N*N];
end
M = inv(M) * N;


% Given D, we need the surface over the center block in the 3x3 grid.
% The matrices U and V are used in the calculation.

u = [N+1:N+N];
u = [u.*u; u; ones(1,N)];
V = u([1 2 1 3 2 1 3 2 3], :)' ;
U = u([1 1 2 1 2 3 2 3 3], :);

z = zigzag(N);
dim=size(X)/N; m=dim(1); n=dim(2);

% Extract the DC coefficients from X.
DC = X(1:N:m*N, 1:N:n*N);

for v = 1:m-2
  vv = v*N;
  for u = 1:n-2
    % Find the coefficients A1--A9.
    A = M * reshape(DC(v:v+2, u:u+2)', 9, 1);
    % Find the surface over the center block of the current 3x3 grid.
    S = dct(V * ( U .* A(:, ones(1, N)) ), N);
    % Now smooth.
    uu = u*N;
    for j = 2:k+1
      y = z(j,1); x = z(j,2);
      if X(vv+y, uu+x) == 0        % this is a candidate for smoothing
	a = S(y,x);                % predicted value
	if abs(a) > Q(y,x)/2       % clamp
	  a =  sign(a) * Q(y,x)/2;
	end
	X(vv+y, uu+x) = a;
      end
    end
  end
end