Image Compression: Reducing Data Redundancies

1. Image Compression Chapter 8

2. 1 Introduction and background The problem: Reducing the amount of data required to represent a digital image. Compression is achieved by removing the data redundancies: Coding redundancy Interpixel redundancy Psychovisual redundancy

3. 1 Introduction and background (cont.) Structral blocks of image compression system. Encoder Decoder The compression ratio where n 1 and n 2 denote the number of information carrying units (usually bits) in the original and encoded images respectively.

4. 1 Introduction and background (cont.)

5. 1 Introduction and background (cont.) function cr = imratio(f1,f2) %IMRATIO Computes the ratio of the bytes in two images / variables. %CR = IMRATIO( F1 , F2 ) returns the ratio of the number of bytes in %variables / files F1 an d F2. If F1 and F2 are an original and compressed %image, respectively, CR is the compression ratio. error(nargchk(2,2,nargin)); %check input arguments cr = bytes(f1) / bytes(f2); %compute the ratio The Function that finds the compression ratio between two images:

6. 1 Introduction and background (cont.) (cont.) % return the number of bytes in input f. If f is a string, assume that it % is an image filename; if not, it is an image variable. function b = bytes(f) if ischar(f) info = dir(f); b = info.bytes; elseif isstruct(f) b = 0; fields = fieldnames(f); for k = 1:length(fields) b = b + bytes(f.(fields{k})); end else info = whos('f'); b = info.bytes; end

7. 1 Introduction and background (cont.) >> r = imratio( (imread( 'bubbles25.jpg' ) ) , b ubbles25.jpg') r = 0.9988 >> f = imread('bubbles.tif'); >> imwrite(f,'bubbles.jpg','jpg') >> r = imratio( (imread( 'bubbles.tif' ) ) , 'bubbles.jpg') r = 14.8578

8. 1 Introduction and background (cont.) Let denote the reconstructed image. Two types of compression Loseless compression: if Loseless compression: if

9. 1 Introduction and background (cont.) In lossy compression the error between is defined by root mean square which is given by

10. 1 Introduction and background (cont.) The M- Function that computes e(rms) and displays both e ( x , y ) and its histogram % COMPARE Computes and displays the error between two matrices. RMSE = COMPARE (F1 , F2, SCALE) returns the root-mean-square error between inputs F1 an d F2, displays a histogram of the difference, and displays a scaled difference image. When SCALE ,s omitted, a scale factor of 1 is used function rmse = compare(f1 , f2 , scale) error(nargchk(2,3,nargin)); if nargin < 3 scale = 1; end %compute the root mean square error e = double(f1) - double(f2); [m,n] = size(e); rmse = sqrt (sum(e(:) .^ 2 ) / (m*n));

11. 1 Introduction and background (cont.) (cont.) %output error image & histogram if an error if rmse %form error histogram. emax= max(abs(e(:))); [h,x]=hist(e(:),emax); if length(h) >= 1 figure; bar(x,h,'k'); %scale teh error image symmetrically and display emax = emax / scale; e = mat2gray(e,[-emax,emax]); figure; imshow(e); end end

12. 1 Introduction and background (cont.) >> r1 = imread('bubbles.tif'); >> r2 = imread('bubbles.jpg'); >> compare(r1, r2,1) >> In E:\matlab\toolbox\images\images\truesize.m (Resize1) at line 302 In E:\matlab\toolbox\images\images\truesize.m at line 40 In E:\matlab\toolbox\images\images\imshow.m at line 168 In E:\matlab\work\matlab_code\compare.m at line 32 ans = 1.5660

13. 1 Introduction and background (cont.) Error histogram

14. 1 Introduction and background (cont.) Error image

15. 2 Coding redundancy Let n k denote the number of times that the k th gray level appears in the image and n denote the total number of pixels in the image. The associated probability for the k th gray level can be expressed as

16. 2 Coding redundancy (cont.) If the number of bits used to represent each value of r k is l ( r k ), then the average number of bits required to represent each pixel is Thus the total number of bits required to code an M N image is MNL avg

17. 2 Coding redundancy (cont.) When fixed variable length that is l ( r k ) = m then

18. 2 Coding redundancy (cont.) The average number of bits requred by Code 2 is

19. 2 Coding redundancy (cont.) How few bits actually are needed to represent the gray levels of an image? Is there a minimum amount of data that is sufficient to describe completely an image without loss information? Information theory provides the mathematical framework to answer these questions.

20. 2 Coding redundancy (cont.) Formulation of generated information Note that if P (E)=1 that is if the event always occurs then I ( E )=0 and no information is attributed to it. No uncertainity associated with the event

21. 2 Coding redundancy (cont.) Given a source of random events from the discrete set of possible events with associated probabilities The average information per source output, called the entropy, is

22. 2 Coding redundancy (cont.) If we assume that the histogram based on gray levels is an estimate of the true probability distribution, the estimate of H can be expressed by

23. 1 Introduction and background (cont.) function h = entropy(x,n) %ENTROPY computes a first-order estimate of the entropy of a matrix. %H = ENTROPY(X,N) returns the first-order estimate of matrix X with N %symbols in bits / symbol. The estimate assumes a statistically independent %source characterized by the relative frequency of occurence of the %elements in X error(nargchk(1,2,nargin)); if nargin < 2 n=256; end

24. 1 Introduction and background (cont.) (cont) x = double(x); xh = hist(x(:),n); xh = xh/sum(xh(:)); % make mask to eliminate 0's since log2(0) = -inf i = find(xh); h = -sum(xh(i) .* log2(xh(i))); %compute entropy

25. 1 Introduction and background (cont.) f = [119 123 168 119;123 119 168 168] ; f = [f; 119 119 107 119 ; 107 107 119 119] ; f = 119 123 168 119 123 119 168 168 119 119 107 119 107 107 119 119 p=hist(f(:),8) p = 3 8 2 0 0 0 0 3

26. 1 Introduction and background (cont.) p = p / sum(p) p = Columns 1 through 7 0.1875 0.5000 0.1250 0 0 0 0 0.1875 H = entropy(f) h = 1.7806

27. Huffman codes Huffman codes are widely used and very effective technique for compressing data. We consider the data to be a sequence of charecters.

28. Huffman codes (cont.) Consider a binary charecter code wherein each charecter is represented by a unique binary string. Fixed-length code: a = 000, b = 001, c = 010, d = 011, e = 100, f = 101 variable-length code: a = 0, b = 101, c = 100, d = 111, e = 1101, f = 1100

29. Huffman codes (cont.) 100 86 14 58 28 14 b:13 a:45 d:16 c:12 f:5 e:9 100 55 25 30 c:12 a:45 a:45 b:13 14 d:16 d:16 f:5 e:9 Fixed-length code 1 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 Variable-length code

30. Huffman codes (cont.) Prefix code: Codes in which no codeword is also a prefix of some other codeword. Encoding for binary code: Example: Variable-length prefix code. a b c Decoding for binary code: Example: Variable-length prefix code.

31. Constructing Huffman codes

32. Constructing Huffman codes Huffmans algorithm assumes that Q is implemented as a binary min-heap. Running time: Line 2 : O ( n ) (uses BUILD-MIN-HEAP ) Line 3-8: O ( n lg n ) (the for loop is executed exactly n- 1 times and each heap operation requires time O ( lg n ) )

33. Constructing Huffman codes: Example c:12 b:13 e:9 f:5 d:16 a:45 14 b:13 c:12 d:16 a:45 f:5 e:9 14 f:5 e:9 d:16 a:45 25 c:12 b:13 0 0 0 0 0 1 1 1

34. Constructing Huffman codes: Example 14 f:5 e:9 d:16 a:45 25 c:12 b:13 30 0 0 0 1 1 1

35. Constructing Huffman codes: Example 14 f:5 e:9 d:16 a:45 25 c:12 b:13 30 55 0 0 0 0 1 1 1 1

36. Constructing Huffman codes: Example 14 f:5 e:9 d:16 a:45 25 c:12 b:13 30 55 100 1 1 1 1 1 0 0 0 0 0

37. Huffman Codes function CODE = huffman(p) %check the input arguments for reasonableness error(nargchk(1,1,nargin)); if (ndims(p) ~= 2) | (min(size(p))>1)| ~isreal(p)|~isnumeric(p) error('P must be a real numeric vector'); end global CODE CODE = cell(length(p),1); if length (p) > 1 p=p /sum(p); s=reduce(p); makecode(s,[]); else CODE ={'1'}; end;

38. Huffman Codes(cont.) %Create a Huffman source reduction tree in a MATLAB cell structure by performing % source symbol reductions until there are only two reduced symbols remaining function s = reduce(p); s= cell(length(p),1) for i=1:length(p) s{i}=i; end while numel(s) > 2 [p,i] = sort(p); p(2) = p(1) + p(2); p(1) = []; s = s(i); s{2} = {s{1},s{2}}; s(1) = []; end

39. Huffman Codes(cont.) %Scan the nodes of a H u ffman source reduction tree recursively to generate %the indicated variable length code words . function makecode(sc,codeword) global CODE if isa(sc,'cell') makecode(sc{1},[codeword 0]); makecode(sc{2},[codeword 1]); else CODE{sc} = char('0'+codeword); end

40. Huffman Codes(cont.) >> p = [0.1875 0.5 0.125 0.1875]; >> c = huffman(p) c = '011' '1' '010' '00'

41. Huffman Encoding >> f2 = uint8 ([2 3 4 2; 3 2 4 4; 2 2 1 2; 1 1 2 2]) f2 = 2 3 4 2 3 2 4 4 2 2 1 2 1 1 2 2 >> whos('f2') Name Size Bytes Class f2 4x4 16 uint8 array Grand total is 16 elements using 16 bytes

42. Huffman Encoding(cont.) >> c = huffman(hist(double(f2(:)),4)) c = '011' '1' '010' '00' >> h1f2 = c(f2(:))' h1f2 = Columns 1 through 9 '1' '010' '1' '011' '010' '1' '1' '011' '00' Columns 10 through 16 '00' '011' '1' '1' '00' '1' '1'

43. Huffman Encoding(cont.) >> whos('h1f2') Name Size Bytes Class h1f2 1x16 1018 cell array Grand total is 45 elements using 1018 bytes >> h2f2 = char (h1f2)' h2f2 = 1010011000011011 1 11 1001 0 0 10 1 1 >> whos('h2f2') Name Size Bytes Class h2f2 3x16 96 char array Grand total is 48 elements using 96 bytes

44. Huffman Encoding(cont.) >> h2f2 = h2f2(:); >> h2f2(h2f2 == ' ') = []; >> whos('h2f2') Name Size Bytes Class h2f2 29x1 58 char array Grand total is 29 elements using 58 bytes

45. Huffman Encoding(cont.) >> h3f2 = mat2huff(f2) h3f2 = size: [4 4] min: 32769 hist: [3 8 2 3] code: [43867 1944] >> whos('h3f2') Name Size Bytes Class h3f2 1x1 518 struct array Grand total is 13 elements using 518 bytes

46. Huffman Encoding(cont.) >> hcode = h3f2.code; >> whos('hcode') Name Size Bytes Class hcode 1x2 4 uint16 array Grand total is 2 elements using 4 bytes >> dec2bin(double(hcode)) ans = 1010101101011011 0000011110011000

47. Huffman Encoding(cont.) >> f = imread('Tracy.tif'); >> c=mat2huff(f); >> cr1 = imratio(f,c) cr1 = 1.2191 >> save SqueezeTracy c; >> cr2 = imratio ('Tracy.tif','SqueezeTracy.mat') cr2 = 1.2627

48. Huffman Decoding function x = huff2mat(y) % HUFF2MAT decodes a Huffman encoded matrix if ~isstruct(y) | ~isfield(y,'min') | ~isfield(y,'size') | ~isfield(y,'hist') | ~isfield(y,'code') error('The input must be a structure as returned by MAT2HUFF'); end sz = double(y.size); m=sz(1); n= sz(2); xmin = double(y.min) - 32768; map = huffman(double(y.hist)); code = cellstr(char('','0','1')); link = [2; 0; 0]; left = [2 3]; found = 0; tofind = length(map);

49. Huffman Encoding(cont.) (cont.) while length(left) & (found < tofind) look = find(strcmp(map, code {left(1)})); if look link(left(1)) = -look; left = left (2:end); found = found + 1; else len = length (code); link(left(1)) = len + 1; link = [link; 0; 0]; code{end + 1} = strcat(code{left(1)},'0'); code{end + 1} = strcat(code{left(1)},'1'); left = left(2:end); left = [left len + 1 len + 2]; end end x = unravel (y.code',link, m*n); %Decode using C 'unravel' x = x + xmin - 1; x = reshape(x,m,n);

50. Huffman Encoding(cont.) (cont.) #include mex.h void unravel (unsigned short *hx, double *link, double *x, double xsz, int hxsz) { int i = 15, j=0, k=0, n=0; while (xsz - k) { if (*(link + n) > 0) { if((*(hx + j) >> i) & 0x0001) n = *(link + n); else n = *(link + n) - 1; if (i) i--; else {j++; i=15;} if( j > hxsz ) mexErrMsgTxt("OutOf code bits??"); }

51. Huffman Encoding(cont.) (cont.) else { *(x + k++) = -*(link + n); n = 0;} } if ( k == xsz - 1 ) *(x + k++) = -*(link + n); } void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { double *link, *x, xsz; unsigned short *hx; int hxsz; if(nrhs != 3) mexErrMsgTxt("Three inputs required.");

52. Huffman Decoding(cont.) (cont.) else if (nlhs > 1) mexErrMsgTxt("Too many output arguments"); if( !mxIsDouble(prhs[2]) || mxIsComplex(prhs[2]) || mxGetN(prhs[2]) * mxGetM(prhs[2]) != 1 ) mexErrMsgTxt("Input XSIZE must be a scalar"); hx = mxGetPr(prhs[0]); link = mxGetPr(prhs[1]); xsz = mxGetScalar(prhs[2]); hxsz = mxGetM(prhs[0]); plhs[0] = mxCreateDoubleMatrix( xsz, 1, mxREAL); x = mxGetPr(plhs[0]); unravel(hx, link, x, xsz, hxsz); }

53. Consider these two images: They have virtually identical histograms The histograms are three modal indicating the presence of three dominant grey levels Interpixel redundancy

54. Interpixel redundancy (cont.) Because the gray levels of the image are not equally probable, variable-length coding can be used to reduce the coding redundancy . Random matches Entropy=7.4253 Compression ratio=1.0704 Aligned matches Entropy=7.3505 Compression ratio=1.0821 Note that the entropy estimates and the compression ratios of the two images are about the same. The variable length coding does not take advantage of the obvious structral relationship between the aligned matches

55. Interpixel redundancy (cont.) In order to reduce pixel redundancies, the 2-D pixel array normally used for human viewing and interpretation must be transformed into a more efficient format. For example the difference between the adjacent pixels can be used to represent an image. Transformations of this type are referred to as mapping . If the original image can be reconstructed from such transformed set is referred to as reversible mappings.

56. Interpixel redundancy (cont.) The following figure shows a lossless predictive coding.

57. Interpixel redundancy (cont.) The following figure shows a lossless predictive coding. The decoder reconstructs the prediction error from the received variable-length code words and performs the inverse operation

58. Interpixel redundancy (cont.) Various methods can be used to generate . In most cases the prediction is formed by a linearcombination of m previous pixels. where m is the order of the predictor and are prediction coefficients . For 1-D linear predictive coding this equation can be rewritten

59. Interpixel redundancy (cont.) Example: Consider encoding the Figure 8.7(c) using the simple linear predictor

60. Psychovisual Redundancy Psychovisual Redundancy is associated with real or quantifiable visual information. Its elimination is desirable because the information itself is not essential for normal visual processing.

61. JPEG compression (cont.) In transform coding, a reversible, linear transform like DFT or DCT is used to map an image into a set of transform coefficients which are then quantized and coded.

62. JPEG compression (cont.) One of the most popular and comprehensive compression standard is the JPEG standard. In the JPEG baseline coding system, which is based on the discrete cosine transform and is adequate for most compression applications. It is performed in four sequential steps:

63. JPEG compression (cont.)

64. JPEG compression (cont.)

65. JPEG compression (cont.) After computing the DCT coefficients, they are normalized in accordance with where , u,v =0,1, ...,7 are the resulting normalized and quantized coefficients, T ( u,v ) DCT of an 8x8 block of image f ( x,y ) and Z ( u,v ) is a transform normalization array like that of Figure 8.12(a)

66. JPEG compression (cont.) After each blocks DCT coefficientsare quantized, the elements of are recorded in accordance with zigzag pattern of Figure 8.12(b). Since the resulting on-dimensionally arrayqualitatively arranged according to increasing spatial frequency,the encoder in Figure 8.11(a) is designed to take advantage of long runs of zero that normally result from the reordering .

67. JPEG 2000 JPEG 2000 is based on the idea that the coefficients of a transform that decorrelates the pixels of an image can be coded more efficiently than the original pixels themselves. If the transform basis functions- wawelets in the JPEG 2000 case- pack most of the important visual information into a small number of coefficients, the remaining coefficients can be quantized coarsely or truncated to zero with little image distortion.

68. JPEG 2000

69. Fractal Image Compression

70. Resim inde ki Benzer Paralar Range : Kk bloklar Domain: Byk bloklar

71. Grntnn bloklara blnmesi Orjinal grnt 4x4 Range Paralar Orjinal 8x8 Domain p aralar 4x4 e indirgenmi halleri

72. Paralarn Benzerlii

73. Paralarn Benzerlii-3 Benzerini arayacamz Range sttedir. Domain ler in orijinal halleri ilk kolonda gsterilmitir. Regresyon modeline gre dzeltilmi domainler ise ikinci stunda gsterilmitir. Yandaki rakamlar residual mean square( rms ) deerlerini gstermektedir. Rms deeri en kk olan domain eldeki range ile eletirilir.

74. Rangeleri Domainler Cinsinden fade Etmek Benzer Range-Domainler bulunduktan sonra bunlar bir dosyada saklanarak resim bilgisini olutururlar. rnek resmimiz 320x200x8 = 512,000 bit (64,000 byte) ieriyor. Domain numaralar 0-999 arasnda olduu iin 10 bitle ifade edilebilir. Parlaklk ve konrast tamsayya evrilip, 0-15 ve 0-31 aralnda kuantize edilirse, 4 ve 5 bitle ifade edilebilir. Dolaysyla, 4000 Range e karlk gelen Domain bilgilerini dosyaya yazarsak : 4000*(Domain Numaras+Parlaklk+Kontrast) = 4000*(10+4+5) = 76,000 bit (9,500 byte) olmaktadr. 512/76 = 6.73:1, (512-76)/512 = %85 sktrmak demektir, stelik %85 sktran dier algoritmalardan ok daha iyi bir kalitede !

75. Grnt Kalitesi D ecompression, siyah bir ekrandan balanarak 8-10 iterasyonla gerekletirilir. Resmin kalitesi sonucun orjinale ne kadar yakn olduuna yani herbir Range-Domain iftinin ne kadar benzer olduk-larna baldr .

76. Decompress ion: Iteration -1

77. Decompress ion: Iteration - 2

78. Decompress ion: Iteration - 3

79. Decompress ion: Iteration -1 0

80. Ayn Resimden Daha Fazla Domain Elde Etmek : Domain Transformasyonlar Rangee uygun Domaini ararken Domainleri olduu gibi brakmayp transformasyon ilemlerine sokabiliriz. Orjinal, 90 o , 180 o ve 270 o derece dndrlr. X eksenine gre yansmalar alnr. Elde edilen 8 ayr Domainin negatifi alnarak, bir Domain-den hareketle 16 farkl Domain elde edilir.