JPEGのアルゴリズムについて

JPEGアルゴリズムの確認

[Reference]
http://en.wikipedia.org/wiki/Jpeg

2006/10/21 平野拓一（東京工業大学）
2006/10/21 Takuichi Hirano (Tokyo Institute of Technology)

BMPファイルの読み込み

SetDirectory["d:/hira2/public_html/hobby/edu/jpeg"]

SetDirectory["d:/Home/takuichi/hira2/public_html/hobby/edu/jpeg"]

g = Import["a.bmp"]

Show[g]

InputForm[g]

ImgList = g[[1, 1]];      

ImgListRed = Map[(#[[1]] &), ImgList, {2}];      
ImgListGreen = Map[(#[[2]] &), ImgList, {2}];      
ImgListBlue = 
            Map[(#[[3]] &), ImgList, {2}];      

nx = 
        Length[ImgList[[1]]];
ny = 
        Length[ImgList];

MakeRGBRasterImage[ImgListRed_, ImgListGreen_, ImgListBlue_] := 
        Module[{}, 
Graphics[
            Raster[Table[{ImgListRed[[j, i]], ImgListGreen[[j, i]], 
                  ImgListBlue[[j, i]]}, {j, 1, ny}, {i, 1, 
                  nx}], 
{{0, 0}, {nx, ny}}, {0, 255}, 
              ColorFunction -> 
                RGBColor]]
];

Show[
      MakeRGBRasterImage[ImgListRed, ImgListGreen, 
        ImgListBlue], 
AspectRatio -> 
        Automatic]

ZeroTable = 
        Table[0, {n, 0, ny - 1}, {m, 0, 
            nx - 1}];
Print["**** Red Component ****"];
Show[
        MakeRGBRasterImage[ImgListRed, ZeroTable, 
          ZeroTable], 
AspectRatio -> 
          Automatic];
Print["**** Green Component ****"];
Show[
        MakeRGBRasterImage[ZeroTable, ImgListGreen, 
          ZeroTable], 
AspectRatio -> 
          Automatic];
Print["**** Blue Component ****"];
Show[
        MakeRGBRasterImage[ZeroTable, ZeroTable, 
          ImgListBlue], 
AspectRatio -> 
          Automatic];

(R,G,B) -> (Y, Cb, Cr)
JPEGではRGBでなく、YCbCr(またはYUVとも言われる)でそれぞれの成分を処理する。YUVはテレビ放送のPAL方式で採用されている色表現。

http://en.wikipedia.org/wiki/YCbCr

The Y component represents the brightness of a pixel 
The Cb and Cr components together represent the chrominance

The human eye can see more detail in the Y component than in Cb and Cr. Using this knowledge, encoders can be designed to compress images more efficiently.

YCbCrFromRGB[{r_, g_, b_}] := {0.299*r + 0.587*g + 
            0.114*b, 
128 - 0.168736*r - 0.331264*g + 
            0.5*b, 
128 + 0.5*r - 0.418688*g - 
            0.081312*b};
RGBFromYCbCr[{y_, cb_, cr_}] := {-179.45579155073017` - 
            1.2175717664097974`*^-6*cb + 1.4019995886573404`*cr + 
            y, 
135.45879483222961` - 
            0.3441356788389385`*cb - 0.7141361555818125`*cr + 
            y, 
-226.81606046360054` + 
            1.7720000660738162`*cb + 4.062980628562637`*^-7*cr + 
            y};

p[0] = {0, 0, 0};
p[1] = {255, 0, 0};
p[2] = {255, 255, 0};
p[3] = {0, 255, 0};
p[4] = {0, 0, 255};
p[5] = {255, 0, 255};
p[6] = {255, 255, 255};
p[7] = {0, 255, 255};
pcubeRGB = {{p[0], p[1]}, {p[1], p[2]}, {p[2], p[3]}, {p[3], 
            p[0]}, 
{p[4], p[5]}, {p[5], p[6]}, {p[6], 
            p[7]}, {p[7], p[4]}, 
{p[0], p[4]}, {p[1], 
            p[5]}, {p[2], p[6]}, {p[3], p[7]}};
pcubeYCbCr = 
        Map[YCbCrFromRGB, pcubeRGB, {2}];
pcubeRGB2 = {{RGBColor[1, 0, 0], 
            Line[{p[0], p[1]}]}, {RGBColor[0, 1, 0], 
            Line[{p[0], p[3]}]}, {RGBColor[0, 0, 1], 
            Line[{p[0], p[4]}]}};
pcubeYCbCr2 = {{RGBColor[1, 0, 0], 
            Line[{YCbCrFromRGB[p[0]], YCbCrFromRGB[p[1]]}]}, {RGBColor[0, 1, 
              0], Line[{YCbCrFromRGB[p[0]], YCbCrFromRGB[p[3]]}]}, {RGBColor[
              0, 0, 1], 
            Line[{YCbCrFromRGB[p[0]], 
                YCbCrFromRGB[p[4]]}]}};
Show[{Graphics3D[{Map[Line, pcubeRGB, 1], 
              pcubeRGB2}], 
Graphics3D[{Map[Line, 
                pcubeYCbCr, 1], 
              pcubeYCbCr2}]}, 
Boxed -> 
          False];

ImgListY = 
        Map[(YCbCrFromRGB[#][[1]] &), ImgList, {2}];      
ImgListCb = 
        Map[(YCbCrFromRGB[#][[2]] &), ImgList, {2}];      
ImgListCr = 
        Map[(YCbCrFromRGB[#][[3]] &), ImgList, {2}];      

ImgListY = ImgListY - 128;

Print["**** Y Component ****"];
ListDensityPlot[
        ImgListY, 
FrameTicks -> 
          False, 
PlotRange -> {-128, 
            127}];
Print["**** Cb Component ****"];
ListDensityPlot[
        ImgListCb, 
FrameTicks -> 
          False, 
PlotRange -> {0, 
            255}];
Print["**** Cr Component ****"];
ListDensityPlot[
        ImgListCr, 
FrameTicks -> 
          False, 
PlotRange -> {0, 
            255}];

DFT (Discrete Fourier Transform)
JPEGではDCTが使われるので、これは関係ないが、基本的なDFTと比較するために計算してる。

FImgListY = Fourier[ImgListY];
FImgListCb = Fourier[ImgListCb];
FImgListCr = Fourier[ImgListCr];

<< Graphics`Legend`;
\(\(\(ShowLegend[
ListDensityPlot[
          Abs[\ Reverse[FImgListY]\ ], 
Mesh -> 
            False, 
PlotRange -> 
            Automatic, 
PlotLabel -> "Amplitude", 
FrameLabel -> {"m, ωx", "n, ωy"}, 
FrameTicks -> 
            None, 
ColorFunction -> \((Hue[\(-0.7\)*\((# - 1)\), 1, 
                  1] &)\), 
DisplayFunction -> 
            Identity]
, {\((Hue[0.7*\((#)\), 1, 1] &)\), 
          16, "Max", "Min", LegendPosition -> {1.1, \(- .4\)}, 
          LegendSize -> {0.4, 1. }, 
          LegendLabel -> ""}
];\)\(
\)
    \)\)
colfun[x_] := 
        RGBColor[0, 3*x, 1] /; (x ≤ 1/3);
colfun[x_] := 
        RGBColor[3*(x - 1/3), 1, 
            1] /; (1/3 < x ≤ 2/3);
\(\(\(colfun[x_] := 
        RGBColor[\(-3\)*\((x - 2/3)\) + 1, \(-3\)*\((x - 2/3)\) + 1, 
            1] /; \((2/3 < x)\);\)\(
\)
    \)\)
\(\(ShowLegend[
ListDensityPlot[
          Arg[Reverse[FImgListY]]*180\/π, 
Mesh -> 
            False, 
PlotRange -> {\(-180\), 
              180}, 
PlotLabel -> "Phase", 
FrameLabel -> {"m, ωx", "n, ωy"}, 
FrameTicks -> 
            None, 
ColorFunction -> 
            colfun, 
DisplayFunction -> 
            Identity]
, {colfun[\(-#\) + 1] &, 
          16, "180", "-180", LegendPosition -> {1.1, \(- .4\)}, 
          LegendSize -> {0.4, 1. }, 
          LegendLabel -> "deg"}
];\)\)

DCT (Discrete Cosine Transform; 離散コサイン変換)

基底の描画（計算と無関係。教育用）

c[p_] := If[p == 0, 1/(Sqrt[2]), 1];
b[u_, v_, m_, n_] := 
        (2/nn)*c[u]*c[v]*Cos[(((2*m + 1)*u*π)/(2*nn))]*
          Cos[(((2*n + 1)*v*π)/(2*nn))];
nn = 8;
buv[u_, v_] := 
        ListDensityPlot[
          Reverse[Table[
              b[u, v, m, n], {m, 0, nn - 1}, {n, 0, 
                nn - 1}]], 
Mesh -> 
            False, 
Frame -> 
            True, 
FrameTicks -> 
            None, 
DisplayFunction -> 
            Identity];
Show[
      GraphicsArray[
        Table[buv[u, v], {u, 0, nn - 1}, {v, 0, 
            nn - 1}]], 
ImageSize -> {320, 
          320}]

２次元リストにDCTを行う関数の定義

[参考]
貴家仁志：よくわかるディジタル画像処理、CQ出版社、p.123

DCT[lis_] := 
        Module[{c, g, nn}, 
c[p_] := 
            If[p == 0, 1/(Sqrt[2]), 1]; 
nn = 
            Length[lis]; 
g[u_, v_] := 
            (2/nn)*c[u]*c[v]*
              Sum[lis[[m + 1, n + 1]]*
                  Cos[(((2*m + 1)*u*π)/(2*nn))]*
                  Cos[(((2*n + 1)*v*π)/(2*nn))], {m, 0, 
                  nn - 1}, {n, 0, nn - 1}]; 
Table[
            g[u, v], {u, 0, nn - 1}, {v, 0, 
              nn - 1}]
];

FImgListY = DCT[ImgListY];
FImgListCb = DCT[ImgListCb];
FImgListCr = DCT[ImgListCr];

    
Print[Round[FImgListY] //  MatrixForm];
Print[
        Round[FImgListCb] //  MatrixForm];
Print[
        Round[FImgListCr] //  MatrixForm];

<< Graphics`Legend`;
\(\(\(ShowLegend[
ListDensityPlot[
          Abs[\ Reverse[FImgListY]\ ], 
Mesh -> 
            False, 
PlotRange -> 
            Automatic, 
PlotLabel -> "Amplitude", 
FrameLabel -> {"m, ωx", "n, ωy"}, 
FrameTicks -> 
            None, 
ColorFunction -> \((Hue[\(-0.7\)*\((# - 1)\), 1, 
                  1] &)\), 
DisplayFunction -> 
            Identity]
, {\((Hue[0.7*\((#)\), 1, 1] &)\), 
          16, "Max", "Min", LegendPosition -> {1.1, \(- .4\)}, 
          LegendSize -> {0.4, 1. }, 
          LegendLabel -> ""}
];\)\(
\)
    \)\)
colfun[x_] := 
        RGBColor[0, 3*x, 1] /; (x ≤ 1/3);
colfun[x_] := 
        RGBColor[3*(x - 1/3), 1, 
            1] /; (1/3 < x ≤ 2/3);
colfun[x_] := 
        RGBColor[-3*(x - 2/3) + 1, -3*(x - 2/3) + 1, 
            1] /; (2/3 < x);

Quantization (量子化)

量子化テーブル

qty=(({{16, 11, 10, 16, 24, 40, 51, 61}, {12, 12, 14, 19, 26, 58, 60, 55}, {14, 13, 16, 24, 40, 57, 69, 56}, {14, 17, 22, 29, 51, 87, 80, 62}, {18, 22, 37, 56, 68, 109, 103, 77}, {24, 35, 55, 64, 81, 104, 113, 92}, {49, 64, 78, 87, 103, 121, 120, 101}, {72, 92, 95, 98, 112, 100, 103, 99}}));      
qtcbcr=(({{17, 18, 24, 47, 99, 99, 99, 99}, {18, 21, 26, 66, 99, 99, 99, 99}, {24, 26, 56, 99, 99, 99, 99, 99}, {47, 66, 99, 99, 99, 99, 99, 99}, {99, 99, 99, 99, 99, 99, 99, 99}, {99, 99, 99, 99, 99, 99, 99, 99}, {99, 99, 99, 99, 99, 99, 99, 99}, {99, 99, 99, 99, 99, 99, 99, 99}}));      ;

text = 
        Table[{Text[
              StyleForm[ToString[Transpose[Reverse[qty]][[i, j]] ], 
                FontSize -> 16, FontWeight -> "Bold", 
                FontColor -> 
                  Hue[0.7*(Transpose[Reverse[qty]][[i, j]]/
                          255. ), 1, 1]], {i, j}]}, {i, 1, 8}, {j, 1, 
            8}];
Show[
        Graphics[{text}], 
PlotRange -> {{0.5, 
              8.5}, {0.5, 8.5}}, 
AspectRatio -> 
          Automatic];

量子化する
sは品質（大きいほど良い）

s = 1;

FImgListY2 = 
        Round[FImgListY/(qty/s)];
FImgListCb2 = 
        Round[FImgListCb/(qtcbcr/s)];
FImgListCr2 = 
        Round[FImgListCr/(qtcbcr/
                s)];

Print[FImgListY2  //  MatrixForm];
Print[FImgListCb2  //  MatrixForm];
Print[FImgListCr2  //  MatrixForm];

本来ならばここでジグザグスキャンを行った後、エントロピー符号化（ハフマン符号化または算術符号化）が入り、JPEGファイルが生成される。

dx = dy = 1. ;
ggrid = 
        Table[{Line[{{dx*(i - 1), 0}, {dx*(i - 1), dy*8}}], 
            Line[{{0, dy*(j - 1)}, {dx*8, dy*(j - 1)}}]}, {i, 1, 
            9}, {j, 1, 9}];
zigzagmat=(({{1, 2, 6, 7, 15, 16, 28, 29}, {3, 5, 8, 14, 17, 27, 30, 43}, {4, 9, 13, 18, 26, 31, 42, 44}, {10, 12, 19, 25, 32, 41, 45, 54}, {11, 20, 24, 33, 40, 46, 53, 55}, {21, 23, 34, 39, 47, 52, 56, 61}, {22, 35, 38, 48, 51, 57, 60, 62}, {36, 37, 49, 50, 58, 59, 63, 64}}));
mat = 
        Transpose[
          Reverse[zigzagmat]];
zigzaglin = 
        Line[Table[Position[mat, i][[1]] - {dx/2, dy/2}, {i, 1, 
              64}]];
text = 
        Table[{Text[
              StyleForm[ToString[mat[[i, j]] ], FontSize -> 16, 
                FontWeight -> "Bold"], {i - dx/2, j - dy/2}]}, {i, 
            1, 8}, {j, 1, 8}];
Show[
        Graphics[{ggrid, {RGBColor[1, 0, 1], AbsoluteThickness[2], 
              zigzaglin}, text}], 
AspectRatio -> 
          Automatic];

ここから先はJPEGファイルの展開過程である。

Inverse Quantization (逆量子化)

FImgListY3 = 
        Round[FImgListY2*(qty/s)];
FImgListCb3 = 
        Round[FImgListCb2*(qtcbcr/s)];
FImgListCr3 = 
        Round[FImgListCr2*(qtcbcr/
                s)];

Print[FImgListY3 //  MatrixForm];
Print[FImgListCb3  //  MatrixForm];
Print[FImgListCr3  //  MatrixForm];

IDCT (Inverse Discrete Cosine Transform; 逆離散コサイン変換)

IDCT[lis_] := 
        Module[{c, f, nn}, 
c[p_] := 
            If[p == 0, 1/(Sqrt[2]), 1]; 
nn = 
            Length[lis]; 
f[m_, n_] := 
            (2/nn)*Sum[
                c[u]*c[v]*lis[[u + 1, v + 1]]*
                  Cos[(((2*m + 1)*u*π)/(2*nn))]*
                  Cos[(((2*n + 1)*v*π)/(2*nn))], {u, 0, 
                  nn - 1}, {v, 0, nn - 1}]; 
Table[
            f[u, v], {u, 0, nn - 1}, {v, 0, 
              nn - 1}]
];

IFFImgListY = IDCT[FImgListY3];
IFFImgListCb = IDCT[FImgListCb3];
IFFImgListCr = IDCT[FImgListCr3];

(Y, Cb, Cr)->(R,G,B)

(Y, Cb, Cr)->(R,G,B)の式の導出（教育用）

Clear[y, cb, cr, r, g, b];
Solve[{y == 
            0.299*r + 0.587*g + 0.114*b, 
cb == 
            128 - 0.168736*r - 0.331264*g + 
              0.5*b, 
cr == 
            128 + 0.5*r - 0.418688*g - 0.081312*b}, {r, g, b}]  //  
      Simplify

(Y, Cb, Cr)->(R,G,B)
を計算

IFFImgListY = IFFImgListY + 128;      
ImgList2 = 
        Table[{IFFImgListY[[i, j]], IFFImgListCb[[i, j]], 
            IFFImgListCr[[i, j]]}, {i, 1, nx}, {j, 1, ny}];

ImgListR2 = 
        Map[(RGBFromYCbCr[#][[1]] &), ImgList2, {2}]  //  
          Round;      
ImgListG2 = 
        Map[(RGBFromYCbCr[#][[2]] &), ImgList2, {2}]  //  
          Round;      
ImgListB2 = 
        Map[(RGBFromYCbCr[#][[3]] &), ImgList2, {2}]  //  
          Round;

画像を表示

Show[
        MakeRGBRasterImage[ImgListR2, ImgListG2, 
          ImgListB2], 
AspectRatio -> 
          Automatic];