diff --git a/src/doc/palette.cpp b/src/doc/palette.cpp
index d609b0c48..35287eab5 100644
--- a/src/doc/palette.cpp
+++ b/src/doc/palette.cpp
@@ -20,11 +20,20 @@
 
 #include <algorithm>
 #include <limits>
+#include <cmath>
 
 namespace doc {
 
 using namespace gfx;
 
+enum class FitCriteria {
+  OLD,
+  RGB,
+  linearizedRGB,
+  CIEXYZ,
+  CIELAB
+};
+
 Palette::Palette()
   : Palette(0, 256)
 {
@@ -321,49 +330,156 @@ void Palette::initBestfit()
   }
 }
 
+// Auxiliary function for rgbToOtherSpace()
+static double f(double t)
+{
+  if (t > 0.00885645171)
+    return std::pow(t, 0.3333333333333333);
+  else
+    return (t / 0.12841855 + 0.137931034);
+}
+
+// Auxiliary function for findBestfit()
+static void rgbToOtherSpace(double& r, double& g, double& b, FitCriteria fc)
+{
+  if (fc == FitCriteria::RGB)
+    return;
+  double Rl, Gl, Bl;
+  // Linearization:
+  r = r / 255.0;
+  g = g / 255.0;
+  b = b / 255.0;
+  if (r <= 0.04045)
+    Rl = r / 12.92;
+  else
+    Rl = std::pow((r + 0.055) / 1.055, 2.4);
+  if (g <= 0.04045)
+    Gl = g / 12.92;
+  else
+    Gl = std::pow((g + 0.055) / 1.055, 2.4);
+  if (b <= 0.04045)
+    Bl = b / 12.92;
+  else
+    Bl = std::pow((b + 0.055) / 1.055, 2.4);
+  if (fc == FitCriteria::linearizedRGB) {
+    r = Rl;
+    g = Gl;
+    b = Bl;
+    return;
+  }
+  // Conversion lineal RGB to CIE XYZ
+  r = 41.24564*Rl + 35.75761 * Gl + 18.04375 * Bl;
+  g = 21.26729*Rl + 71.51522 * Gl + 7.2175   * Bl;
+  b = 1.93339*Rl  + 11.91920 * Gl + 95.03041 * Bl;
+  switch (fc) {
+
+    case FitCriteria::CIEXYZ:
+      return;
+
+    case FitCriteria::CIELAB: {
+      // Converting CIEXYZ to CIELAB:
+      // For Standard Illuminant D65:
+      //  const double xn = 95.0489;
+      //  const double yn = 100.0;
+      //  const double zn = 108.884;
+      double xxn = r / 95.0489;
+      double yyn = g / 100.0;
+      double zzn = b / 108.884;
+      double fyyn = f(yyn);
+
+      double Lstar = 116.0 * fyyn - 16.0;
+      double aStar = 500.0 * (f(xxn) - fyyn);
+      double bStar = 200.0 * (fyyn - f(zzn));
+
+      r = Lstar;
+      g = aStar;
+      b = bStar;
+      return;
+    }
+  }
+}
+
 int Palette::findBestfit(int r, int g, int b, int a, int mask_index) const
 {
   ASSERT(r >= 0 && r <= 255);
   ASSERT(g >= 0 && g <= 255);
   ASSERT(b >= 0 && b <= 255);
   ASSERT(a >= 0 && a <= 255);
-  ASSERT(!col_diff.empty());
 
-  r >>= 3;
-  g >>= 3;
-  b >>= 3;
-  a >>= 3;
+  FitCriteria fc = FitCriteria::OLD;
 
-  // Mask index is like alpha = 0, so we can use it as transparent color.
-  if (a == 0 && mask_index >= 0)
-    return mask_index;
+  if (fc == FitCriteria::OLD) {
+    ASSERT(!col_diff.empty());
 
-  int bestfit = 0;
-  int lowest = std::numeric_limits<int>::max();
-  int size = std::min(256, int(m_colors.size()));
+    r >>= 3;
+    g >>= 3;
+    b >>= 3;
+    a >>= 3;
 
-  for (int i=0; i<size; ++i) {
-    color_t rgb = m_colors[i];
+    // Mask index is like alpha = 0, so we can use it as transparent color.
+    if (a == 0 && mask_index >= 0)
+      return mask_index;
 
-    int coldiff = col_diff_g[((rgba_getg(rgb)>>3) - g) & 127];
-    if (coldiff < lowest) {
-      coldiff += col_diff_r[(((rgba_getr(rgb)>>3) - r) & 127)];
+    int bestfit = 0;
+    int lowest = std::numeric_limits<int>::max();
+    int size = std::min(256, int(m_colors.size()));
+
+    for (int i=0; i<size; ++i) {
+      color_t rgb = m_colors[i];
+
+      int coldiff = col_diff_g[((rgba_getg(rgb)>>3) - g) & 127];
       if (coldiff < lowest) {
-        coldiff += col_diff_b[(((rgba_getb(rgb)>>3) - b) & 127)];
+        coldiff += col_diff_r[(((rgba_getr(rgb)>>3) - r) & 127)];
         if (coldiff < lowest) {
-          coldiff += col_diff_a[(((rgba_geta(rgb)>>3) - a) & 127)];
-          if (coldiff < lowest && i != mask_index) {
-            if (coldiff == 0)
-              return i;
+          coldiff += col_diff_b[(((rgba_getb(rgb)>>3) - b) & 127)];
+          if (coldiff < lowest) {
+            coldiff += col_diff_a[(((rgba_geta(rgb)>>3) - a) & 127)];
+            if (coldiff < lowest && i != mask_index) {
+              if (coldiff == 0)
+                return i;
 
-            bestfit = i;
-            lowest = coldiff;
+              bestfit = i;
+              lowest = coldiff;
+            }
           }
         }
       }
     }
+
+    return bestfit;
   }
 
+  if (a == 0 && mask_index >= 0)
+    return mask_index;
+
+  int bestfit = 0;
+  double lowest = std::numeric_limits<double>::max();
+  int size = m_colors.size();
+  // Linearice:
+  double x = double(r);
+  double y = double(g);
+  double z = double(b);
+
+  rgbToOtherSpace(x, y, z, fc);
+
+  for (int i=0; i<size; ++i) {
+    color_t rgb = m_colors[i];
+    double Xpal = double(rgba_getr(rgb));
+    double Ypal = double(rgba_getg(rgb));
+    double Zpal = double(rgba_getb(rgb));
+    // Palette color conversion RGB-->XYZ and r,g,b is assumed CIE XYZ
+    rgbToOtherSpace(Xpal, Ypal, Zpal, fc);
+    double xDiff = x - Xpal;
+    double yDiff = y - Ypal;
+    double zDiff = z - Zpal;
+    double aDiff = double(a - rgba_geta(rgb)) / 128.0;
+
+    double diff = xDiff * xDiff + yDiff * yDiff + zDiff * zDiff + aDiff * aDiff;
+    if (diff < lowest) {
+      lowest = diff;
+      bestfit = i;
+    }
+  }
   return bestfit;
 }