dumb by default


dir: Home / Posts / Polyfitting Colormap
published-date: 25 May 2025 21:16 +0700
categories: [quickie] [python-hijinks]
tags: [python]

Polyfitting Colormap


This one’s a quick writeup on how to encode/decode matplotlib colormap (CC0 license for due diligence) while avoiding using external dependencies.

Why not use a library or pasting predefined lookup table you may ask? With polyfit we can reduce lookup table into only n-ths degree polynomial coefficients! And generate lookup table (or even sample singular value) at compile time or at boot! Also I want to show other uses case of polyfit outside statistics.

You can obtain cmap lookup table from the following link. The following example uses inferno variant.

https://github.com/BIDS/colormap/blob/master/colormaps.py

  1import numpy as np
  2import matplotlib.pyplot as plt
  3
  4N_POLY = 4
  5data = np.array(_inferno_data)
  6xs = np.linspace(0, 1, 512)
  7
  8# generator function to decode sample at array x from polynomial coef array
  9def fp(x, p):
 10  def fp_i(v, p):
 11    s = 0
 12    n = len(p)
 13    for i in range(n):
 14      s += (v ** (n - i - 1)) * p[i]
 15    return s
 16  buf = [0 for _ in range(len(x))]
 17  for i in range(len(x)):
 18    buf[i] = fp_i(x[i], p)
 19  return buf
 20
 21# cleaner numpy version
 22def fp_np(x, p):
 23  n = np.flip(np.arange(len(p)))
 24  a = np.expand_dims(np.array(x), axis=1)
 25  a = a ** n
 26  a = a * p
 27  a = a.sum(axis=1)
 28  a[a<0] = 0
 29  a[a>1] = 1
 30  return a
 31
 32xx = np.linspace(0,1,data.shape[0])
 33p_r = np.polyfit(xx, data[:, 0], N_POLY)
 34p_g = np.polyfit(xx, data[:, 1], N_POLY)
 35p_b = np.polyfit(xx, data[:, 2], N_POLY)
 36
 37repr_ps = lambda p: '[' + ', '.join(map(lambda v: f'{v:.16e}', p)) + ']'
 38str_p_r = repr_ps(p_r)
 39str_p_g = repr_ps(p_g)
 40str_p_b = repr_ps(p_b)
 41
 42print('copy the following each channel p coef:')
 43print('p_r =', str_p_r)
 44print('p_g =', str_p_g)
 45print('p_b =', str_p_b)
 46
 47# the rest is optional, but allows us to see how decoded
 48# values may differ from source.
 49
 50# here we re-parse from string we generate previously 
 51# to simulate possible precision losses from pasted values.
 52p_r = eval(str_p_r)
 53p_g = eval(str_p_g)
 54p_b = eval(str_p_b)
 55
 56src = np.vstack([
 57  np.interp(xs, xx, data[:, 0]),
 58  np.interp(xs, xx, data[:, 1]),
 59  np.interp(xs, xx, data[:, 2])
 60]).T
 61
 62dec = np.vstack([
 63  fp_np(xs, p_r),
 64  fp_np(xs, p_g),
 65  fp_np(xs, p_b)
 66]).T
 67
 68diff = src - dec
 69diff01 = (diff + 1) / 2
 70err = np.sqrt(np.mean(diff**2))
 71
 72fig = plt.figure(figsize=(6, 5))
 73
 74sfig = fig.subfigures(3, 1, height_ratios=[2, 2, 1])
 75
 76xt = np.linspace(0, 1, 32)
 77for i in range(3):
 78  ax = sfig[0].add_subplot(1, 3, i+1)
 79  ax.scatter(xt, fp_np(xt, [p_r, p_g, p_b][i]), c='black', marker='.')
 80  ax.plot(xs, src[:, i], c=['r', 'g', 'b'][i], linewidth=2)
 81  if i>0: ax.get_yaxis().set_visible(False)
 82  ax.set_ylim((-.01, 1.01))
 83  ax.set_title(['red', 'green', 'blue'][i])
 84
 85ax = sfig[1].add_subplot(3, 1, 1)
 86ax.imshow(np.vstack([[src], [src]]), aspect='auto')
 87ax.text(-.02, .5, 'source', va='center', ha='right', transform=ax.transAxes)
 88ax.set_axis_off()
 89
 90ax = sfig[1].add_subplot(3, 1, 3)
 91ax.imshow(np.vstack([[dec], [dec]]), aspect='auto')
 92ax.text(-.02, .5, 'poly\ndecoded', va='center', ha='right', transform=ax.transAxes)
 93ax.set_axis_off()
 94
 95ax = sfig[1].add_subplot(3, 1, 2)
 96ax.imshow(np.vstack([[diff01], [diff01]]), aspect='auto')
 97ax.text(-.02, .5, 'diff', va='center', ha='right', transform=ax.transAxes)
 98ax.set_axis_off()
 99
100ax = sfig[2].add_subplot()
101ax.set_title(f'diff each channel, rmse:{err:.02e}')
102ax.plot([0, 1], [0, 0], c='black', linestyle=':')
103ax.plot([0, 1], [err, err], c='black', linestyle=':')
104ax.plot([0, 1], [-err, -err], c='black', linestyle=':')
105ax.plot(xs, diff[:, 0], c='r')
106ax.plot(xs, diff[:, 1], c='g')
107ax.plot(xs, diff[:, 2], c='b')
108ax.set_xlim((xs[0], xs[-1]))
109ax.grid()
110
111fig.show()

The code above will create this graph.

image

The first three graph shows each channel value distribution with the dotted lines being the decoded values. Here we can see polyfit does not “fit” fully to sampled source.

The second three are the rendered version with diff. If you look closely with only 4 coef you can see bands of lightblue tint.

The third one visualizes each channel differences.

With inferno color table, which contains large amount of values, can be represented with only these coef (assuming we generate only 4).

p_r = [1.6557762312002666e+00, -5.0765628176037012e+00, 3.6114437522600484e+00, 7.5501780705212251e-01, -1.0159038301211065e-02]
p_g = [-2.0627850606099467e+00, 4.4300104345314564e+00, -1.8488324386347577e+00, 5.0539609607001779e-01, -1.5872300948496232e-03]
p_b = [7.0917066940883133e+00, -6.6777021038402200e+00, -2.3104640769293701e+00, 2.5883795816309192e+00, 2.6268954622100660e-03]

You can minimize the error increasing the number of coef globally or on each channel on line 33-35

Tangent: you can create or edit matplotlib colormap using viscm.





Built with Hugo | previoip (c) 2025