Color Image Quality Assessment Based on CIEDE2000, Yang Yang, Jun Ming and Nenghai Yu, 2012

analyze_ciede2000.py

#!/usr/bin/env python3.4
from collections import deque
import sys
import numpy as np
from skimage import color
import random
from scipy import ndimage

# Assuming BT.709
yuv2rgb = np.array([
    [1., 0., 1.28033], [1., -0.21482, -0.38059], [1., 2.12798, 0.]
])


def generate_lab_from_ycbcr():
    Y_i = range(0, 257, 4)
    C_i = range(0, 257, 4)
    Y_buf = [x for x in Y_i for y in C_i for z in C_i]
    Cb_buf = [y for x in Y_i for y in C_i for z in C_i]
    Cr_buf = [z for x in Y_i for y in C_i for z in C_i]
    W, H = len(Y_buf), 1
    Y = (np.array(Y_buf).reshape((H, W)) - 16.) / 219.
    Cb = (np.array(Cb_buf).reshape((H, W)) - 128.) / 224.
    Cr = (np.array(Cr_buf).reshape((H, W)) - 128.) / 224.
    YCbCr = np.dstack((Y, Cb, Cr))
    rgb = np.dot(YCbCr, yuv2rgb.T)
    lab = color.rgb2lab(rgb) * 128
    L = lab[:, :, 0].reshape((65, 65, 65))
    a = lab[:, :, 1].reshape((65, 65, 65))
    b = lab[:, :, 2].reshape((65, 65, 65))
    return np.array((L, a, b), dtype=np.int16)


lab_from_ycbcr = generate_lab_from_ycbcr()


def ycbcr2lab(img):
    coords = img.reshape((img.shape[0] * img.shape[1], 3)).T * (1. / 4.)
    L = ndimage.map_coordinates(lab_from_ycbcr[0], coords, order=1)
    a = ndimage.map_coordinates(lab_from_ycbcr[1], coords, order=1)
    b = ndimage.map_coordinates(lab_from_ycbcr[2], coords, order=1)
    return np.dstack((L, a, b)).reshape(img.shape) * (1. / 128.)


img = np.array(
    [[x, y, z] for x in range(256) for y in range(256) for z in range(256)],
    dtype=np.uint8).reshape((4096, 4096, 3))
inter = ycbcr2lab(img).flatten()
ref = (color.rgb2lab(np.dot(
    [[((x - 16) / 219., (y - 128) / 224., (z - 128) / 224.)
      for x in range(256) for y in range(256) for z in range(256)]
     ], yuv2rgb.T))).flatten()
print(ref)
print(inter)
print(np.square((inter - ref)).mean())

visualize_ciede2000.py

#!/usr/bin/env python3.4
import math
import sys
import y4m
import numpy as np
from collections import deque
from skimage import color
import matplotlib.pyplot as plt
from matplotlib import mlab, cm
from scipy import ndimage

# Assuming BT.709

yuv2rgb = np.array([
  [1.,  0.,       1.28033],
  [1., -0.21482, -0.38059],
  [1.,  2.12798,  0.]])
box2 = np.ones((2,2))

def decode_y4m_buffer(frame):
  width = frame.headers['W']
  height = frame.headers['H']
  Y = np.ndarray(shape=(height,width), buffer=frame.buffer, dtype='uint8')
  Cb = np.ndarray(shape=(height//2,width//2), buffer=frame.buffer, offset=width*height, dtype='uint8')
  Cr = np.ndarray(shape=(height//2,width//2), buffer=frame.buffer, offset=int(width*height*1.25), dtype='uint8')
  return np.dot(np.dstack(((Y - 16.) / 219., np.kron((Cb - 128.) / 224., box2), np.kron((Cr - 128.) / 224., box2))), yuv2rgb.T)


# See "Color Image Quality Assessment Based on CIEDE2000" Yang Yang, Jun Ming and Nenghai Yu, 2012
# http://dx.doi.org/10.1155/2012/273723

def process_pair(ref, recons, recons2):
  ref_img = decode_y4m_buffer(ref)
  recons_img = decode_y4m_buffer(recons)
  recons2_img = decode_y4m_buffer(recons2)
  ref_lab = color.rgb2lab(ref_img)
  recons_lab = color.rgb2lab(recons_img)
  recons2_lab = color.rgb2lab(recons2_img)
  dE = color.deltaE_ciede2000(ref_lab, recons_lab, kL=0.65, kC=1.0, kH=4.0)
  dE2 = color.deltaE_ciede2000(ref_lab, recons2_lab, kL=0.65, kC=1.0, kH=4.0)
  peak = 48.
  enhanced = ndimage.gaussian_filter(np.maximum(dE, dE2) - dE2, sigma=5)
  dE = np.clip(dE, 0., peak)/peak
  fig = plt.figure(figsize=(ref.headers['W']/96., ref.headers['H']/96.), dpi=96, frameon=False)
  ax=fig.add_subplot(1,1,1)
  plt.axis('off')
  plt.imshow(np.power(dE, 2./3.), cmap=cm.Greys_r)
  levels = np.linspace(enhanced.mean(), enhanced.max(), 10)
  plt.contourf(enhanced, levels=levels, cmap=cm.jet, alpha=0.3, linewidth=0.)
  fig.subplots_adjust(bottom = 0)
  fig.subplots_adjust(top = 1)
  fig.subplots_adjust(right = 1)
  fig.subplots_adjust(left = 0)
  extent = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
  fig.savefig('%03d.png' % ref.count, bbox_inches=extent, dpi=96)
  plt.close()
  print('%03d.png' % ref.count)
  pass

ref_frames = deque()
recons_frames = deque()
recons2_frames = deque()

def process_ref(frame):
  ref_frames.append(frame)
  if recons_frames and recons2_frames:
    process_pair(ref_frames.popleft(), recons_frames.popleft(), recons2_frames.popleft())
  pass

def process_recons(frame):
  recons_frames.append(frame)
  if ref_frames and recons2_frames:
    process_pair(ref_frames.popleft(), recons_frames.popleft(), recons2_frames.popleft())
  pass

def process_recons2(frame):
  recons2_frames.append(frame)
  if ref_frames and recons_frames:
    process_pair(ref_frames.popleft(), recons_frames.popleft(), recons2_frames.popleft())
  pass

def main(args):
  ref_parser = y4m.Reader(process_ref)
  recons_parser = y4m.Reader(process_recons)
  recons_parser2 = y4m.Reader(process_recons2)
  with open(args[1], 'r') as ref:
    with open(args[2], 'r') as recons:
      with open(args[3], 'r') as recons2:
        while True:
          data = ref.buffer.read(2*1024*1024)
          if not data: break
          ref_parser.decode(data)
          data = recons.buffer.read(2*1024*1024)
          if not data: break
          recons_parser.decode(data)
          data = recons2.buffer.read(2*1024*1024)
          if not data: break
          recons_parser2.decode(data)

if __name__ == '__main__':
  main(sys.argv)

y4m_ciede2000.py

#!/usr/bin/env python3
import math
import sys
import y4m
import numpy as np
from collections import deque
from skimage import color

# Assuming BT.709

yuv2rgb = np.array([
  [1.,  0.,       1.28033],
  [1., -0.21482, -0.38059],
  [1.,  2.12798,  0.]])
box2 = np.ones((2,2))

def decode_y4m_buffer(frame):
  width = frame.headers['W']
  height = frame.headers['H']
  Y = np.ndarray(shape=(height,width), buffer=frame.buffer, dtype='uint8')
  Cb = np.ndarray(shape=(height//2,width//2), buffer=frame.buffer, offset=width*height, dtype='uint8')
  Cr = np.ndarray(shape=(height//2,width//2), buffer=frame.buffer, offset=int(width*height*1.25), dtype='uint8')
  return np.dot(np.dstack(((Y - 16.) / 219., np.kron((Cb - 128.) / 224., box2), np.kron((Cr - 128.) / 224., box2))), yuv2rgb.T)


# See "Color Image Quality Assessment Based on CIEDE2000" Yang Yang, Jun Ming and Nenghai Yu, 2012
# http://dx.doi.org/10.1155/2012/273723

scores = []

def process_pair(ref, recons):
  ref_img = decode_y4m_buffer(ref)
  recons_img = decode_y4m_buffer(recons)
  ref_lab = color.rgb2lab(ref_img)
  recons_lab = color.rgb2lab(recons_img)
  dE = color.deltaE_ciede2000(ref_lab, recons_lab, kL=0.65, kC=1.0, kH=4.0).mean()
  scores.append(math.log(24./dE, 2.))
  print('%03d: %f' % ( ref.count, scores[-1] ))
  pass

ref_frames = deque()
recons_frames = deque()

def process_ref(frame):
  ref_frames.append(frame)
  if recons_frames:
    process_pair(ref_frames.popleft(), recons_frames.popleft())
  pass

def process_recons(frame):
  recons_frames.append(frame)
  if ref_frames:
    process_pair(ref_frames.popleft(), recons_frames.popleft())
  pass

def main(args):
  ref_parser = y4m.Reader(process_ref)
  recons_parser = y4m.Reader(process_recons)
  with open(args[1], 'r') as ref:
    with open(args[2], 'r') as recons:
      while True:
        data = ref.buffer.read(2*1024*1024)
        if not data: break
        ref_parser.decode(data)
        data = recons.buffer.read(2*1024*1024)
        if not data: break
        recons_parser.decode(data)
  print('AVG: %f' % np.array(scores).mean())
  

if __name__ == '__main__':
  main(sys.argv)