LINUX.ORG.RU
решено Форум — Development  

2D Discrete Cosine Transform + inverse 2D DCT


0

1

Предупреждаю - в теории этого преобразования я практически ничего не понимаю, поэтому вопрос может показаться тупым. Пробую выполнить FDCT над матрицей 8x8 со значениями 255. http://en.wikipedia.org/wiki/Discrete_cosine_transform#DCT-II Затем - выполнить над результатом inverse DCT. Но оригинальная 8x8 матрица при этом не восстанавливается. Как выполнить обратное преобразование посредством IDCT? Вот программа на Python: 

#!/usr/bin/python

# pixels
MorigB = [
[255,255,255,255,255,255,255,255],
[255,255,255,255,255,255,255,255],
[255,255,255,255,255,255,255,255],
[255,255,255,255,255,255,255,255],
[255,255,255,255,255,255,255,255],
[255,255,255,255,255,255,255,255],
[255,255,255,255,255,255,255,255],
[255,255,255,255,255,255,255,255],
]

from math import sqrt,cos,pi
from sys import stdout
M = 8
N = 8

def dct(a):
  F = []
  for u in xrange(N):
    F = F + [[]]
    for v in xrange(M):
      s = 0.0
      lu = 1.0
      if (u == 0): lu = 1/float(sqrt(2))
      lv = 1.0
      if (v == 0): lv = 1/float(sqrt(2))
      for i in xrange(N):
        for j in xrange(M):
          d = (a[i][j]);
          d *= cos(((pi*u)*(2*i+1))/float(2*N))
          d *= cos(((pi*v)*(2*j+1))/float(2*M))
          s += d
      s *= lu
      s *= lv
      s *= sqrt(2.0/N)
      s *= sqrt(2.0/M)
      F[u] = F[u] + [s,]
  return F

def idct(a):
  F = []
  for u in xrange(N):
    F = F + [[]]
    for v in xrange(M):
      s = 0.0
      for i in xrange(N):
        for j in xrange(M):
          li = 1
          if (i == 0): li = 1/float(sqrt(2))
          lj = 1
          if (j == 0): lj = 1/float(sqrt(2))
          d = (a[i][j]);
          d *= cos(((pi*u)*(2*i+1))/float(2*N))
          d *= cos(((pi*v)*(2*j+1))/float(2*M))
          d *= (li*lj)
          s += d
      s *= (1/4.0)
      F[u] = F[u] + [s,]
  return F

def add(a,d):
  F = []
  for u in xrange(N):
    F = F + [[]]
    for v in xrange(M):
      F[u] = F[u] + [a[u][v]+d,]
  return F

def print_(a):
  print "---"
  for u in xrange(N):
    for v in xrange(M):
      stdout.write(" %6.0f" % (int(a[u][v])))
    stdout.write("\n")

a = add(MorigB,-128)
print_(a)
b = dct(a)
print_(b)
c = idct(b)
d = add(c,+128)
print_(d)


$ ./dct.py 
---
    127    127    127    127    127    127    127    127
    127    127    127    127    127    127    127    127
    127    127    127    127    127    127    127    127
    127    127    127    127    127    127    127    127
    127    127    127    127    127    127    127    127
    127    127    127    127    127    127    127    127
    127    127    127    127    127    127    127    127
    127    127    127    127    127    127    127    127
---
   1015      0      0      0      0      0      0      0
      0      0      0      0      0      0      0      0
      0      0      0      0      0      0      0      0
      0      0      0      0      0      0      0      0
      0      0      0      0      0      0      0      0
      0      0      0      0      0      0      0      0
      0      0      0      0      0      0      0      0
      0      0      0      0      0      0      0      0
---
    255    252    245    233    217    198    176    152
    252    250    243    231    216    197    175    152
    245    243    236    225    210    193    172    150
    233    231    225    215    202    186    168    148
    217    216    210    202    191    177    162    145
    198    197    193    186    177    167    155    141
    176    175    172    168    162    155    146    137
    152    152    150    148    145    141    137    132

★★★★★

Последнее исправление: pacify (всего исправлений: 1)
← Растровое изображение по полярным координатам точек в OpenGL
[Tcl/Tk][Serial][Нипонял] Еще вопрос по fileevent →

Ну кто ж так пишет на питоне.

По теме - ты бы потренировался на кош^Wодномерных данных, для начала. Раз не понимаешь, что к чему...

anonymous
()
Ответ на: комментарий от fang

Ну, он не настолько lossy, чтобы так сильно результаты отличались от первоначальных. Сдается мне, программная реализация хромает (что-нибудь где-нибудь неправильно округляется, float'ы используются и т.п.).

А вообще, конечно, чтобы было полное совпадение с оригиналом, надо или Фурье или вейвлеты использовать...

Eddy_Em ☆☆☆☆☆
()

idct:

li = 1
if (i == 0): li = 1/float(sqrt(2))
lj = 1
if (j == 0): lj = 1/float(sqrt(2))

Здесь надо умножать на обратный коэффициент, чтобы вернуть постоянную составляющую в первоначальный вид.

s *= (1/4.0)

Откуда ты это взял? Надо так же, как и при прямом преобразовании, умножать на sqrt(2.0/N) * sqrt(2.0/M), чтобы выполнялось условие равенства нормирующего множителя при двустороннем преобразовании (2/N) * (2/M).

fang ★
()
Ответ на: комментарий от Eddy_Em

> А вообще, конечно, чтобы было полное совпадение с оригиналом

Кэп подсказывает, что, поскольку преобразование обратимо, то полнота совпадения с оригиналом зависит только от точности представления чисел.

или вейвлеты использовать


Я давно заметил: для тебя вейвлеты — буквально навязчивая идея какая-то. Всюду их пихаешь, куда надо и не надо.

fang ★
()
Ответ на: комментарий от fang

Кэп подсказывает, что, поскольку преобразование обратимо, то полнота совпадения с оригиналом зависит только от точности представления чисел.

В принципе - да. Лучше float не использовать здесь. Я уже как-то приводил пример «чудес», возникающих при работе с float (когда [A*B-A] не одно и то же, что [A*(B-1.f)]).

Я давно заметил: для тебя вейвлеты — буквально навязчивая идея какая-то. Всюду их пихаешь, куда надо и не надо.

Ну нравятся они мне - что поделать-то? =)

Тебе, может, тоже какие-то преобразования нравятся, которые ты пихаешь куда надо, и куда не надо? :D

Eddy_Em ☆☆☆☆☆
()
Ответ на: комментарий от fang

Здесь надо умножать на обратный коэффициент, чтобы вернуть постоянную составляющую в первоначальный вид.

Так, что-ли?

def idct(a):
  F = []
  for u in xrange(N):
    F = F + [[]]
    for v in xrange(M):
      s = 0.0
      for i in xrange(N):
        for j in xrange(M):
          li = 1
          if (i == 0): li = float(sqrt(2))
          lj = 1
          if (j == 0): lj = float(sqrt(2))
          d = (a[i][j]);
          d *= cos(((pi*u)*(2*i+1))/float(2*N))
          d *= cos(((pi*v)*(2*j+1))/float(2*M))
          d *= (li*lj)
          s += d
      s *= sqrt(2.0/N)
      s *= sqrt(2.0/M)
      F[u] = F[u] + [s,]
  return F

Все-равно, матрица неравномерная получается - сомножитель-то действует только на DC-коэффициент (0,0). 

a = add(MorigB,-128)
b = dct(a)
print_(idct(b))
    508    498    469    422    359    282    194     99
    498    488    460    414    352    276    190     97
    469    460    433    390    331    260    179     91
    422    414    390    351    298    234    161     82
    359    352    331    298    253    199    137     70
    282    276    260    234    199    156    108     55
    194    190    179    161    137    108     74     37
     99     97     91     82     70     55     37     19

s *= (1/4.0) Откуда ты это взял? Надо ... умножать на sqrt(2.0/N) * sqrt(2.0/M)

 — а это и равно 1/4 ... (просто я эти формулы скопипастил из разных источников, чаще всего там упрощения всякие типа 1/4).

pacify ★★★★★
() автор топика
Ответ на: комментарий от fang

И поэтому в первом же элементе вместо 127 получается 255.

127 = 255 - 128 см. первую строчку, где из исходной матрицы вычитается 128: 

a = add(MorigB,-128)
print_(a)
b = dct(a)
print_(b)
c = idct(b)
d = add(c,+128)
print_(d)
Т.е. dct/idct должно восстановить равномерную матрицу из +127 (а потом я добавляю снова +128).

P.S. И все-равно, я не понимаю - как сумма косинусов может восстановить константу? Даже приближенно.

pacify ★★★★★
() автор топика
Ответ на: комментарий от fang

Кэп подсказывает, что, поскольку преобразование обратимо,
то полнота совпадения с оригиналом зависит только от точности представления чисел.

Я попробовал заменить в idct() пару строчек: 

>          d *= cos(((pi*u)*(2*i+1))/float(2*N))
>          d *= cos(((pi*v)*(2*j+1))/float(2*M))

          d *= cos(((pi*u)*(2*i))/float(2*N))
          d *= cos(((pi*v)*(2*j))/float(2*M))
(что подсказала логика - для константной функции болжен учитываться только (0,0) DCT-коэффициент);

Результат: 
---
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
---
 1015.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00
   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00
   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00
   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00
   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00
   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00
   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00
   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00
---
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
 126.00 127.00 126.00 127.00 127.00 126.00 127.00 126.00
 126.00 127.00 126.00 127.00 127.00 126.00 127.00 126.00
 126.00 127.00 126.00 127.00 127.00 126.00 127.00 126.00
 126.00 126.00 126.00 126.00 126.00 126.00 126.00 126.00
 126.00 126.00 126.00 126.00 126.00 126.00 126.00 126.00
 126.00 126.00 126.00 126.00 126.00 126.00 126.00 126.00

Но на других данных: 

---
 -128.00 -128.00 -128.00 127.00 127.00 127.00 127.00 127.00
 -128.00 -128.00 -128.00 127.00 127.00 127.00 127.00 127.00
 -128.00 -128.00 -128.00 127.00 127.00 127.00 127.00 127.00
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
---
 729.00 -320.00 -124.00  46.00  95.00  31.00 -51.00 -63.00
 -320.00 -357.00 -139.00  51.00 106.00  34.00 -57.00 -71.00
 -124.00 -139.00 -54.00  20.00  41.00  13.00 -22.00 -27.00
  46.00  51.00  20.00  -7.00 -15.00  -5.00   8.00  10.00
  95.00 106.00  41.00 -15.00 -31.00 -10.00  17.00  21.00
  31.00  34.00  13.00  -5.00 -10.00  -3.00   5.00   6.00
 -51.00 -57.00 -22.00   8.00  17.00   5.00  -9.00 -11.00
 -63.00 -71.00 -27.00  10.00  21.00   6.00 -11.00 -14.00
---
 -176.00 -122.00 -194.00  -8.00 162.00 107.00 141.00 114.00
 -122.00 -78.00 -137.00  15.00 156.00 110.00 138.00 117.00
 -194.00 -137.00 -213.00 -16.00 164.00 105.00 142.00 114.00
  -8.00  15.00 -16.00  66.00 142.00 118.00 133.00 121.00
 162.00 156.00 164.00 142.00 122.00 129.00 125.00 128.00
 107.00 110.00 105.00 118.00 129.00 125.00 127.00 126.00
 141.00 138.00 142.00 133.00 125.00 127.00 126.00 127.00
 114.00 117.00 114.00 121.00 128.00 126.00 127.00 126.00
Результат неудовлетворительный.

Сейчас буду пробовать idct в Матлабе.

pacify ★★★★★
() автор топика
Ответ на: комментарий от pacify

Какой заказ ? у них уже все готово

The WebM video hardware IP designs support encoding and decoding WebM/VP8 video up to 1080p resolution. Implemented as RTL source code (VHDL or Verilog), both WebM IP designs are available under a no-cost license.

Deliverables

RTL source code (VHDL or Verilog) with RTL test bench and test data.
ANSI C source code for hardware drivers, with software test bench and test data.
Technical documentation: hardware and software integration guides and application programming interface (API) manuals.
Reference C-model (bit-accurate with RTL).

bukinist
()
Ответ на: комментарий от pacify

Главная ошибка — в перепутанной индексации

d *= cos(((pi*u)*(2*i+1))/float(2*N))
d *= cos(((pi*v)*(2*j+1))/float(2*M))

Вместо этого должно быть

d *= cos(((pi*i)*(2*u+1))/float(2*N))
d *= cos(((pi*j)*(2*v+1))/float(2*M))

Либо в условиях циклов сделай замену u <-> i, v <-> j.

fang ★
()
Ответ на: комментарий от fang

Главная ошибка — в перепутанной индексации


Вау! Спасибо. Всё заработало.

Какой заказ ? у них уже все готово


Нас не все устраивает, из того, что они предлагают.

P.S. Всё, вопрос закрыт, всем спасибо.

Пример: 

---
 -128.00 -128.00 -128.00 127.00 127.00 127.00 127.00 127.00
 -128.00 -128.00 -128.00 127.00 127.00 127.00 127.00 127.00
 -128.00 -128.00 -128.00 127.00 127.00 127.00 127.00 127.00
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
 127.00 127.00 127.00 127.00 127.00 127.00 127.00 127.00
---
 729.00 -320.00 -124.00  46.00  95.00  31.00 -51.00 -63.00
 -320.00 -357.00 -139.00  51.00 106.00  34.00 -57.00 -71.00
 -124.00 -139.00 -54.00  20.00  41.00  13.00 -22.00 -27.00
  46.00  51.00  20.00  -7.00 -15.00  -5.00   8.00  10.00
  95.00 106.00  41.00 -15.00 -31.00 -10.00  17.00  21.00
  31.00  34.00  13.00  -5.00 -10.00  -3.00   5.00   6.00
 -51.00 -57.00 -22.00   8.00  17.00   5.00  -9.00 -11.00
 -63.00 -71.00 -27.00  10.00  21.00   6.00 -11.00 -14.00
---
 -128.00 -127.00 -128.00 127.00 126.00 127.00 126.00 126.00
 -128.00 -128.00 -128.00 127.00 126.00 127.00 127.00 126.00
 -128.00 -128.00 -128.00 127.00 126.00 126.00 126.00 126.00
 127.00 127.00 127.00 127.00 127.00 126.00 127.00 126.00
 126.00 126.00 126.00 127.00 126.00 126.00 126.00 126.00
 126.00 127.00 127.00 127.00 126.00 126.00 126.00 126.00
 126.00 126.00 127.00 127.00 126.00 127.00 127.00 126.00
 126.00 126.00 126.00 126.00 126.00 126.00 126.00 126.00

pacify ★★★★★
() автор топика
Ответ на: комментарий от anonymous

> Ну кто ж так пишет на питоне.

А что не так с кодом?
Стандарты какие-то не соблюдаю?

Просто моими первыми языками были: ассемблер для ZX, BASIC (ZX) и Turbo Pascal.

Поэтому привычней писать простенькие программки на Питоне. Руби (когда я выбирал скриптовый язык) был еще непопулярен.

pacify ★★★★★
() автор топика
Ответ на: комментарий от pacify

Ты используешь решения «влоб» заместо использования идиоматичных решений на питоне.

Например, твоя функция add() записывается как «def add(c,x): return map(lambda x: map(lambda x:x+c, x), x)». Без двух вложенных циклов, причем циклов через xrange().

Ты используешь часто такой подход: «for i in xrange(N): <process> x[i]», даже когда тебе индекс не нужен и можно обойтись «for i in x».

Да и чего только стоит твоя матрица 8х8 (ты прям мастер паттерна копи-паст). Её можно записать так: [[255]*8]*8

ky-san
()
Ответ на: комментарий от ky-san

Просто мне надо больше наглядности, чтобы показать другим, не понимающим в программировании особо.

А за [[255]*8]*8 - спасибо, буду использовать. Я такое только для строк использовал: " «*level (для отступов).

pacify ★★★★★
() автор топика
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.
← Растровое изображение по полярным координатам точек в OpenGL
Development
[Tcl/Tk][Serial][Нипонял] Еще вопрос по fileevent →