Стандарты сжатия изображения с потерей качества

Страницы работы

Содержание работы

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ 
ИНФОРМАЦИОННЫХ ТЕХНОЛОГИЙ, МЕХАНИКИ И ОПТИКИ

Отчет по практикуму
«Стандарты сжатия изображения с потерей качества.»

задание 3

Студенты  Чиликин.А.И

Дыбленко А.В

Харитонов С.А

Преподаватель  Бочарова И. Е.

Санкт-Петербург

2009

Текст задания

  1. Выполнить преобразование формата изображения RGB в YUV.
  2. Выполнить ДКП для компоненты Y цветного изображения. Проквантовать полученные коэффициенты. Выполнить преобразование коэффициентов постоянного тока в последовательность категорий и амплитуд.
  3. Исследовать статистическое характеристики полученного потока данных и оценить возможный коэффициент сжатия изображения.
  4. Оценить величину отношения сигнал-шум. Получить восстановленное изображение и вывести его на экран.
  5. Сопоставить полученные результаты с характеристиками стандартного JPEG-коедра.

Исходное изображение:

Листингкода MatLab

function main()

  path = 'd:\files\study\media\3\';

  SourceBinFile = ReadFile([path,'girl.bmp']);

  % заголовки и дата

  bfOffBits = SourceBinFile(11:14);

  biWidth = SourceBinFile(19:22);

  biHeight = SourceBinFile(23:26);

  DATA = SourceBinFile(bfOffBits+1:numel(SourceBinFile));

  headers = SourceBinFile(1:bfOffBits);

  % цвета

  Np = biWidth(1)*biHeight(1);%ширину*высоту

  bpos = [1:3:Np*3];

  BData = DATA(bpos);

  gpos = [2:3:Np*3];

  GData = DATA(gpos);

  rpos = [3:3:Np*3];

  RData = DATA(rpos);

  % RED

  img = [headers];

  Newimg(rpos) = RData;

  img = [img; Newimg'];

  WriteFile(path, 'rad.bmp', img);

  img = [];

  Newimg = [];

  % GREEN

  img = [headers];

  Newimg(gpos) = GData;

  img = [img; Newimg'];

  img = [img; 0];

  WriteFile(path, 'green.bmp', img);

  img = [];

  Newimg = [];

  %  BLUE

  img = [headers];

  Newimg(bpos) = BData;

  img = [img; Newimg'];

  img = [img; 0; 0];

  WriteFile(path, 'blue.bmp', img);

  img = [];

  Newimg = [];

  % RGB 2 YUV

  YData = round(0.299*RData + 0.578*GData + 0.114*BData);

  UData = round((BData - YData)*0.5643 + 128);

  VData = round((RData - YData)*0.7132 + 128);

  % Y-составляющая

  GYUVData = YData;

  RYUVData = YData;

  BYUVData = YData;

  img = [headers];

  Newimg(rpos) = RYUVData;

  Newimg(gpos) = GYUVData;

  Newimg(bpos) = BYUVData;

  img = [img; Newimg'];

  WriteFile(path, 'y.bmp', img);

  img = [];

  Newimg = [];

  % U-составляющая

  GYUVData = UData;

  RYUVData = UData;

  BYUVData = UData;

  img = [headers];

  Newimg(rpos) = RYUVData;

  Newimg(gpos) = GYUVData;

  Newimg(bpos) = BYUVData;

  img = [img; Newimg'];

  WriteFile(path, 'u.bmp', img);

  img = [];

  Newimg = [];

  % V-составляющая

  GYUVData = VData;

  RYUVData = VData;

  BYUVData = VData;

  img = [headers];

  Newimg(rpos) = RYUVData;

  Newimg(gpos) = GYUVData;

  Newimg(bpos) = BYUVData;

  img = [img; Newimg'];

  WriteFile(path, 'v.bmp', img);

  img = [];

  Newimg = [];

  % ДКП

  % блоки 8х8 пикселей

  width = biWidth(1)/8;

  height = biHeight(1)/8;

  % КПЕТ

  AC = [];

  % КПОТ

  DC = [];

  for row = 1:height

    for column = 1:width

      % блок

      Block(1, :) = YData( ((row-1)*width*64 + (column-1)*8 +1):((row-1)*width*64 + (column)*8) )';

      Block(2, :) = YData( ((row-1)*width*64 + width*8 + (column-1)*8 +1):((row-1)*width*64  + width*8 + (column)*8) )';

      Block(3, :) = YData( ((row-1)*width*64 + 2*width*8 + (column-1)*8 +1):((row-1)*width*64 + 2*width*8 + (column)*8) )';

      Block(4, :) = YData( ((row-1)*width*64 + 3*width*8 + (column-1)*8 +1):((row-1)*width*64 + 3*width*8 + (column)*8) )';

      Block(5, :) = YData( ((row-1)*width*64 + 4*width*8 + (column-1)*8 +1):((row-1)*width*64 + 4*width*8 + (column)*8) )';

      Block(6, :) = YData( ((row-1)*width*64 + 5*width*8 + (column-1)*8 +1):((row-1)*width*64 + 5*width*8 + (column)*8) )';

      Block(7, :) = YData( ((row-1)*width*64 + 6*width*8 + (column-1)*8 +1):((row-1)*width*64 + 6*width*8 + (column)*8) )';

      Block(8, :) = YData( ((row-1)*width*64 + 7*width*8 + (column-1)*8 +1):((row-1)*width*64 + 7*width*8 + (column)*8) )';

      DCT = Dct(Block);

      Q = [10,12,9,15,24,40,51,61

          10,11,13,19,26,58,60,55

          15,16,16,24,40,57,69,56

          12,17,22,29,51,87,80,62

          18,22,37,56,68,100,103,77

          24,35,55,64,81,104,112,92

          49,64,78,87,103,121,120,101

          72,92,95,98,112,100,103,99];

      % Квантование

      for i=1:8

        for j=1:8

          QDCT(i, j) = fix(DCT(i, j)/Q(i, j));

        end

      end

      % Сохраним коэффициенты постоянного тока

      if((row == 1)&&(column == 1))

        DC(1) = QDCT(1, 1);

      else

        DC((row-1)*width + column) = QDCT(1, 1) - DC((row-1)*width + column - 1);

      end

      % коэффициенты переменного тока z

      Zway = [2, 9, 17, 10, 3, 4, 11, 18, 25, 33, 26, 19, 12, 5, 6, 13, 20, 27, 34, 41, 49, 42, 35, 28, 21, 14, 7, 8, 15, 22, 29, 36, 43, 50, 57, 58, 51, 44, 37, 30, 23, 16, 24, 31, 38, 45, 52, 59, 60, 53, 46, 39, 32, 40, 47, 54, 61, 62, 55, 48, 56, 63, 64];

      temp = QDCT';

      ACthread = temp(Zway);

      % потоки серий и длины

      i = 1;

      l = 0;

      curAC = [];

      map((row-1)*width + column) = 1;

      while(i < 64)

        if(ACthread(i) == 0)

          l = l+1;

        else

          map((row-1)*width + column) = 0;

          curAC = [curAC, [l; ACthread(i)]];

        end

        i = i + 1;

      end

      % запишем в общий поток с учетом карты блоков

      if(map((row-1)*width + column) == 0 )

        AC = [AC, curAC];

      end;

      % восстановим квантованный блок

      for i=1:8

        for j=1:8

          RDCT(i, j) = QDCT(i, j)*Q(i, j);

        end

      end

      RBlock = UnDct(RDCT);

      RYData( ((row-1)*width*64 + (column-1)*8 +1):((row-1)*width*64 + (column)*8) ) = RBlock(1, :)';

      RYData( ((row-1)*width*64 + width*8 + (column-1)*8 +1):((row-1)*width*64  + width*8 + (column)*8) ) = RBlock(2, :)';

      RYData( ((row-1)*width*64 + 2*width*8 + (column-1)*8 +1):((row-1)*width*64 + 2*width*8 + (column)*8) ) = RBlock(3, :)';

      RYData( ((row-1)*width*64 + 3*width*8 + (column-1)*8 +1):((row-1)*width*64 + 3*width*8 + (column)*8) ) = RBlock(4, :)';

      RYData( ((row-1)*width*64 + 4*width*8 + (column-1)*8 +1):((row-1)*width*64 + 4*width*8 + (column)*8) ) = RBlock(5, :)';

      RYData( ((row-1)*width*64 + 5*width*8 + (column-1)*8 +1):((row-1)*width*64 + 5*width*8 + (column)*8) ) = RBlock(6, :)';

      RYData( ((row-1)*width*64 + 6*width*8 + (column-1)*8 +1):((row-1)*width*64 + 6*width*8 + (column)*8) ) = RBlock(7, :)';

      RYData( ((row-1)*width*64 + 7*width*8 + (column-1)*8 +1):((row-1)*width*64 + 7*width*8 + (column)*8) ) = RBlock(8, :)'; 

    end

  end

  %количество байт для хранения DC

  for i=1:numel(DC)

    %считаем количество байт на каждую амплитуду

    ADC(i) = size(dec2bin(abs(DC(i))), 2);

  end

  % Энтропия

  for i=min(ADC):max(ADC)

    probD(i) = sum(ADC == i)/numel(ADC);

  end

  eDC = 0;

  for i=min(ADC):max(ADC)

    eDC = eDC + probD(i)*log2(probD(i));

  end

  eDC = -eDC;

  disp('----------start-------------');

  disp(['Entropy = ', num2str(eDC)]);

  % Кол-во бит на всю последовательность DC

  bDC = sum(ADC) + eDC*numel(ADC);

  disp(['Bit count = ', num2str(bDC)]);

  % Посчитаем сколько бит надо, чтобы закодировать последовательность AC

  % для начала высчитаем для амплитуд

  for i=1:numel(AC(2, :))

    % считаем количество байт на каждую амплитуду

    AAC(i) = size(dec2bin(abs(AC(2, i))), 2);

  end

  for i=min(AAC):max(AAC)

    probA(i) = sum(AAC == i)/numel(AAC);

  end

  eAAC = 0;

  for i=min(AAC):max(AAC)

    if(probA(i) ~= 0) 

      eAAC = eAAC + probA(i)*log2(probA(i));

    end

  end

  eAAC = -eAAC;

  % теперь для длин серий

  for i=min(AC(1, :)):max(AC(1, :))

    probAC(i+1) = sum(AC(1, :) == i)/numel(AC(1, :));

  end

  eAC = 0;

  for i=min(AC(1, :)):max(AC(1, :))

    if(probAC(i+1) ~= 0)

      eAC = eAC + probAC(i+1)*log2(probAC(i+1));

    end

  end

  eAC = -eAC;

  % количество бит на хранение AC

  bAC = sum(AAC) + numel(AC(1, :))*eAC + numel(AC(2, :))*eAAC;

  disp(['AC sequence size = ', num2str(bAC)]);

  % учтем величину бит на хранение карты ненулевых блоков

  bmap = numel(map);

  disp(['Bit count for non-null block = ', num2str(bmap)]);

  % оценим сжатие

  Koeff = 8*width*height*64/(bAC + bDC + bmap);

  disp(['Compression = ', num2str(Koeff)]);

  % вычислим отношение сигнал/шум

  % шум

  Ena = 0;

  for i=1:numel(YData)

    Ena = Ena + (YData(i) - RYData(i))*(YData(i) - RYData(i));

  end

  Ena = Ena/(height*width*64);

  PSNR = 10*log10((255*255)/Ena);

  disp(['Signal/Noise: ', num2str(PSNR)]);

  disp('-----------end--------------');

  RYData = RYData';

  %YUV 2 RGB

  GYUVData = round(RYData - 0.714*(VData - 128) - 0.334*(UData - 128));

  RYUVData = round(RYData + 1.402*(VData - 128));

  BYUVData = round(RYData + 1.772*(UData - 128));

  img = [headers];

  Newimg(rpos) = RYUVData;

  Newimg(gpos) = GYUVData;

  Newimg(bpos) = BYUVData;

  img = [img; Newimg'];

  WriteFile(path, 'recovery.bmp', img);

  img = [];

  Newimg = [];

end

Полученные результаты

Энтропия: 2.5263

Количество Бит: 2386.8856

Размер всей последовательности переменных токов: 13352.1274

Величину бит на хранение карты ненулевых блоков: 300

Оценка Сжатия: 9.5766 (при сжатие по jpeg c качеством 90% составляет 6.5)

Сигнал/Шум: 31.1904

Восстановленное изображение:

Разбитые на цвета:

Y-составляющая:

U-составляющая:

V-составляющая:

Похожие материалы

Информация о работе