deblockc

jpeg的matlab代码

# deblockc - calculate AC coefficients in a smoothing scheme
#
# SYNOPSIS
#	deblockc()
#	deblockc(k)
#	deblockc(k, N)
#
# DESCRIPTION
#       deblockc is a Maple procedure to find the AC coefficients (in terms 
#       of the DC coefficients) which could be used as part of a smoothing 
#       procedure to reduce blocking artifacts (the telltale discontinuities
#       between blocks which often follow aggressive quantizing) in JPEG.
#
#       The polynomial
#          p = a_1 x^2 y^2 + a_2 x^2 y + a_3 x y^2 + a_4 x^2 + a_5 x y
#                + a_6 y^2 + a_7 x +a_8 y + a_9
#       is fitted using the DC coefficients from nearest neighbors:
#
#           ------------------- x
#          | DC1 | DC2 | DC3 |
#          |-----------------|
#          | DC4 | DC5 | DC6 |
#          |-----------------|
#          | DC7 | DC8 | DC9 |
#          |-----------------
#          y
#
#       The coefficients are chosen so that the average over the NxN 
#       DC_i subblock is DC_i. (There is a scaling factor--see below.) 
#       N defaults to 8; k specifies the number of coefficients to 
#       find (defaults to 5).
#
#       The idea is that low-frequency zero AC coefficients in the center 
#       block will be will be replaced by values from the (transformed)
#       polynomial. The procedure is described in Pennebaker and Mitchell,
#       "JPEG Still Image Data Compression Standard", Van Nostrand Reinhold.
#
# EXAMPLE
#       The two statements
#         AC := deblockc():
#         AC := deblockc(5,8): 
#       produce the same results. The zigzag sequence AC[0,0], AC[0,1],
#       AC[1,0], ..., AC[0,3] is displayed (first symbolically, then
#       numerically), and the results are stored in AC.
#
# BUGS
#       This is a hack.

# 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.

deblockc := proc()
  local N, terms, i, X, AC, row, col, v,u, y, x, n, z,
        p, dct, zigzag;

  if nargs < 2 then N := 8; else N := args[2] fi;
  if nargs < 1 then terms := 5; else terms := args[1] fi;

for i from 1 to 9 do; 
  a.i := evaln(a.i);
od;

p := proc(y,x)
 a1*x^2*y^2 + a2*x^2*y + a3*x*y^2 + a4*x^2 + a5*x*y + a6*y^2 
 + a7*x +a8*y + a9
end:

# In the "standard" transform, the DC coefficient is N^2/(2*sqrt (2N))
# times the block averages. Adjustment made below.
for i from 1 to 9 do;
  row := floor((i-1)/3) * N; col := ((i-1) mod 3) * N;
  eq.i := sum(sum(p(y,x),y=row..row+N-1), 
                 x=col..col+N-1) - DC.i:
od:

assign( solve({seq(eq.i, i=1..9)}, {seq(a.i, i=1..9)}) );

X := array(0..N-1, 0..N-1);
for y from 0 to N-1 do;
  for x from 0 to N-1 do;
     X[y,x] := p(y+N, x+N)
  od
od;

dct := proc(X, v, u)
  local x, y, w, N, C;
  N := op(2,eval(X))[1]; N := op(2,N) - op(1,N) + 1;
  C := proc(w)
    if not type(w, 'numeric') then
       RETURN( 'procname(args)' )
    fi;
    if w = 0 then 1/sqrt(2) else 1 fi
  end:
  w := Pi/(2*N);
  C(v)*C(u)/sqrt(2*N) * 
  sum(cos((2*y+1)*v*w) * sum(X[y,x]*cos((2*x+1)*u*w), x=0..N-1), 
      y=0..N-1)
end:

zigzag := proc(N)
  local u, v, inc, n, z;
  z := array(0..N*N-1, 0..1);
  u := -1; v := 1; inc := 1;

  for n from 0 to N*N-1 do; 
    v := v - inc; u := u + inc;
    if (u >= N) then
      v := v + 2; u := N-1; inc := -1;
    elif (v = -1) then
      v := 0; inc := -1;
    elif (v >= N) then
      u := u + 2; v := N-1; inc := 1;
    elif (u = -1) then 
      u := 0; inc := 1;
    fi;
    z[n, 0] := v; z[n,1] := u;
  od;
  z
end:

z := zigzag(N); AC := array(0..N, 0..N);

for n from 0 to terms do;
  v := z[n,0]; u := z[n,1];
  AC[v,u] :=  collect( dct(X,v,u), {seq(DC.i, i=1..9)} ) * 2*sqrt(2*N);
  print(evaln(AC[v,u]) = AC[v,u]);
od:

for n from 0 to terms do;
  v := z[n,0]; u := z[n,1];
  print(evaln(AC[v,u]) = evalf(AC[v,u]));
od:

AC
end:

`help/text/deblockc` := TEXT(
`FUNCTION: deblockc - find coefficients in deblocking scheme`,
`   `,
`CALLING SEQUENCE:`,
`   deblock(k, N)`,
`   `,
`PARAMETERS:`,
`   k  - (optional, default 5) number of terms to find`,
`   N  - (optional, default 8) size of block in the transform`,
`   `,
`SYNOPSIS:`,
`- The function deblockc finds coefficients of a polynomial to be used`,
`  by 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 matrix of coefficients is returned.`,
` `,
`- See Pennebaker and Mitchell, "JPEG Still Image Data Compression Standard"`,
`  Van Nostrand Reinhold, 1993.`
):

#save `deblockc.m`;
#quit