#!/usr/bin/env python # coding: utf-8
# In[1]: get_ipython().run_line_magic('matplotlib', 'inline') import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import PchipInterpolator raw_data = np.loadtxt('out.csv', delimiter=',') ms = sorted(set(raw_data.T[0])) qs = sorted(set(raw_data.T[1])) M_interp = 16 Ms = np.linspace(min(ms), max(ms), (len(ms) - 1)*M_interp + 1) def aggregate(data, drop=None, label=None): pixels = data.T[5]*data.T[6]*3//2 for Q in qs: rate = [] psnr = [] de = [] for M in ms: idx = (data.T[0] == M) & (data.T[1] == Q) if drop is not None and label is not None: idx = idx & (label != drop) px_idx = pixels[idx] px_sum = px_idx.sum() rate.append(data[idx,2].sum()*8/px_sum) psnr.append((data[idx,3] * px_idx).sum()/px_sum) de.append((data[idx,4] * px_idx).sum()/px_sum) # Interpolate the results rate_ = np.poly1d(np.polyfit(ms, np.log(rate), 3)) psnr_ = np.poly1d(np.polyfit(ms, psnr, 3)) de_ = np.poly1d(np.polyfit(ms, de, 3)) for M in Ms: yield (M, Q, np.exp(rate_(M)), psnr_(M), de_(M)) data = np.array([v for v in aggregate(raw_data)]) # # Quality vs bitrate (color vs brightness) # # ### Analysis by [David Michael Barr](https://github.com/barrbrain) # # [CIEDE2000](https://en.wikipedia.org/wiki/Color_difference#CIEDE2000) is a perceptually uniform measure of color difference. # # [PSNR](https://en.wikipedia.org/wiki/Peak_signal-to-noise_ratio) is an essential measure of error, in this case restricted to differences in brightness. # # In the following plots, we compare the trends in quality vs bitrate for various allocations of color detail. # In[2]: def plot_rd(D, N, C): X = data.T[2][data.T[0] == N] Y = data.T[D][data.T[0] == N] Y_ = np.linspace(Y.min(), Y.max(), 1000) X_ = np.exp(PchipInterpolator(Y[::-1], np.log(X[::-1]))(Y_)) plt.semilogx(X_, Y_, color=C, lw=.5) plt.figure(figsize=(16, 8)) plt.subplot(121) plt.title('subset3') plt.xlabel('bits per pixel') plt.ylabel('CIEDE2000') for N in range(40, 101): plot_rd(4, N, plt.cm.gray(N*3-100)) plot_rd(4, 68, 'r') plt.xlim([data.T[2].min(), data.T[2].max()]) plt.ylim([data.T[4].min(), data.T[4].max()]); plt.subplot(122) plt.title('subset3') plt.xlabel('bits per pixel') plt.ylabel('PSNR') for N in range(40, 101): plot_rd(3, N, plt.cm.gray(N*3-100)) plot_rd(3, 68, 'r') plt.xlim([data.T[2].min(), data.T[2].max()]) plt.ylim([data.T[3].min(), data.T[3].max()]); # We found the balance of color detail that maximized perceived color matching but that ended up significantly harming details in brightness for small improvements in color accuracy. So we traded the marginal gains in color accuracy for greater losses in brightness detail by stepping back to the point where there is an even trade between the two. # In[3]: def get_bdr(D): top = min([data.T[D][data.T[0] == N].max() for N in Ms]) bottom = max([data.T[D][data.T[0] == N].min() for N in Ms]) def interp(N): Y = np.log(data.T[2][data.T[0] == N]) X = data.T[D][data.T[0] == N] pchip = PchipInterpolator(X[::-1], Y[::-1], extrapolate=True) return pchip.integrate(bottom, top) / (top - bottom) base = interp(96) return [np.exp(interp(N) - base) - 1 for N in Ms] de_bdr = get_bdr(4) de_bdr_ = np.polyfit(x=Ms, y=de_bdr, deg=3) de_bdr__ = np.polyder(de_bdr_) psnr_bdr = get_bdr(3) psnr_bdr_ = np.polyfit(x=Ms, y=psnr_bdr, deg=3) psnr_bdr__ = np.polyder(psnr_bdr_) best = np.roots(de_bdr__+psnr_bdr__)[0] plt.figure(figsize=(16, 9)) plt.title('subset3') plt.xlabel('chroma quantizer gradient Q6') plt.ylabel('BD-rate') # Q8 to Q6 (/4), since we sampled with a step of 4. plt.scatter(Ms[::M_interp]/4, de_bdr[::M_interp], s=3) plt.scatter(Ms[::M_interp]/4, psnr_bdr[::M_interp], s=3) plt.plot(Ms/4, np.poly1d(de_bdr_)(Ms), label='CIEDE2000') plt.plot(Ms/4, np.poly1d(psnr_bdr_)(Ms), label='PSNR') plt.vlines(x=best/4, ymax=max(de_bdr), ymin=min(psnr_bdr), label='best - '+str(int(round(best/4)))) plt.xlim([10, 25]) plt.ylim([min(psnr_bdr), max(de_bdr)]) plt.legend(); np.poly1d(de_bdr_)(best), np.poly1d(psnr_bdr_)(best), best*16 # We repeated the analysis many times with slightly different collections of images to show the typical balance point and how much it varied according to which images were tested. # In[4]: best = [] raw_label = raw_data.T[7].astype(int) - 1 label_map = np.arange(raw_label.max() + 1) % 4 import multiprocessing def optimize(drop): data = np.array([v for v in aggregate(raw_data, drop, label)]) def get_bdr(D): top = min([data.T[D][data.T[0] == N].max() for N in Ms]) bottom = max([data.T[D][data.T[0] == N].min() for N in Ms]) def interp(N): Y = np.log(data.T[2][data.T[0] == N]) X = data.T[D][data.T[0] == N] pchip = PchipInterpolator(X[::-1], Y[::-1], extrapolate=True) return pchip.integrate(bottom, top) / (top - bottom) base = interp(96) return [np.exp(interp(N) - base) - 1 for N in Ms] de_bdr = get_bdr(4) de_bdr_ = np.polyfit(x=Ms, y=de_bdr, deg=3) de_bdr__ = np.polyder(de_bdr_) psnr_bdr = get_bdr(3) psnr_bdr_ = np.polyfit(x=Ms, y=psnr_bdr, deg=3) psnr_bdr__ = np.polyder(psnr_bdr_) return np.roots(de_bdr__+psnr_bdr__)[0] for drop_outter in range(25): np.random.shuffle(label_map) label = label_map[raw_label] pool = multiprocessing.Pool(4) best.extend(pool.map(optimize, range(4))) pool.close() pool.join() # Q8 to Q6 (/4), since we sampled with a step of 4. best = np.array(sorted(set(best)))/4 best # In[5]: import seaborn as sns from scipy.stats import mstats from scipy.stats import norm best_ = mstats.winsorize(best,(.05,.05)) plt.figure(figsize=(16, 9)) sns.distplot(best_, rug=True, hist=False, fit=norm, kde=True); best_.mean(), best_.std()