如何进一步优化此色差功能？答案

【问题标题】：How can I further optimize this color difference function?如何进一步优化此色差功能？
【发布时间】：2011-03-04 05:09:54
【问题描述】：

我已经制作了这个函数来计算 CIE Lab 色彩空间中的色差，但它缺乏速度。由于我不是Java专家，想知道周围有没有Java高手有一些技巧可以提高这里的速度。

代码基于注释块中提到的matlab函数。

/**
 * Compute the CIEDE2000 color-difference between the sample color with
 * CIELab coordinates 'sample' and a standard color with CIELab coordinates
 * 'std'
 *
 * Based on the article:
 * "The CIEDE2000 Color-Difference Formula: Implementation Notes,
 * Supplementary Test Data, and Mathematical Observations,", G. Sharma,
 * W. Wu, E. N. Dalal, submitted to Color Research and Application,
 * January 2004.
 * available at http://www.ece.rochester.edu/~gsharma/ciede2000/
 */
public static double deltaE2000(double[] lab1, double[] lab2)
{
    double L1 = lab1[0];
    double a1 = lab1[1];
    double b1 = lab1[2];

    double L2 = lab2[0];
    double a2 = lab2[1];
    double b2 = lab2[2];

    // Cab = sqrt(a^2 + b^2)
    double Cab1 = Math.sqrt(a1 * a1 + b1 * b1);
    double Cab2 = Math.sqrt(a2 * a2 + b2 * b2);

    // CabAvg = (Cab1 + Cab2) / 2
    double CabAvg = (Cab1 + Cab2) / 2;

    // G = 1 + (1 - sqrt((CabAvg^7) / (CabAvg^7 + 25^7))) / 2
    double CabAvg7 = Math.pow(CabAvg, 7);
    double G = 1 + (1 - Math.sqrt(CabAvg7 / (CabAvg7 + 6103515625.0))) / 2;

    // ap = G * a
    double ap1 = G * a1;
    double ap2 = G * a2;

    // Cp = sqrt(ap^2 + b^2)
    double Cp1 = Math.sqrt(ap1 * ap1 + b1 * b1);
    double Cp2 = Math.sqrt(ap2 * ap2 + b2 * b2);

    // CpProd = (Cp1 * Cp2)
    double CpProd = Cp1 * Cp2;

    // hp1 = atan2(b1, ap1)
    double hp1 = Math.atan2(b1, ap1);
    // ensure hue is between 0 and 2pi
    if (hp1 < 0) {
        // hp1 = hp1 + 2pi
        hp1 += 6.283185307179586476925286766559;
    }

    // hp2 = atan2(b2, ap2)
    double hp2 = Math.atan2(b2, ap2);
    // ensure hue is between 0 and 2pi
    if (hp2 < 0) {
        // hp2 = hp2 + 2pi
        hp2 += 6.283185307179586476925286766559;
    }

    // dL = L2 - L1
    double dL = L2 - L1;

    // dC = Cp2 - Cp1
    double dC = Cp2 - Cp1;

    // computation of hue difference
    double dhp = 0.0;
    // set hue difference to zero if the product of chromas is zero
    if (CpProd != 0) {
        // dhp = hp2 - hp1
        dhp = hp2 - hp1;
        if (dhp > Math.PI) {
            // dhp = dhp - 2pi
            dhp -= 6.283185307179586476925286766559;
        } else if (dhp < -Math.PI) {
            // dhp = dhp + 2pi
            dhp += 6.283185307179586476925286766559;
        }
    }

    // dH = 2 * sqrt(CpProd) * sin(dhp / 2)
    double dH = 2 * Math.sqrt(CpProd) * Math.sin(dhp / 2);

    // weighting functions
    // Lp = (L1 + L2) / 2 - 50
    double Lp = (L1 + L2) / 2 - 50;

    // Cp = (Cp1 + Cp2) / 2
    double Cp = (Cp1 + Cp2) / 2;

    // average hue computation
    // hp = (hp1 + hp2) / 2
    double hp = (hp1 + hp2) / 2;

    // identify positions for which abs hue diff exceeds 180 degrees
    if (Math.abs(hp1 - hp2) > Math.PI) {
        // hp = hp - pi
        hp -= Math.PI;
    }
    // ensure hue is between 0 and 2pi
    if (hp < 0) {
        // hp = hp + 2pi
        hp += 6.283185307179586476925286766559;
    }

    // LpSqr = Lp^2
    double LpSqr = Lp * Lp;

    // Sl = 1 + 0.015 * LpSqr / sqrt(20 + LpSqr)
    double Sl = 1 + 0.015 * LpSqr / Math.sqrt(20 + LpSqr);

    // Sc = 1 + 0.045 * Cp
    double Sc = 1 + 0.045 * Cp;

    // T = 1 - 0.17 * cos(hp - pi / 6) +
    //       + 0.24 * cos(2 * hp) +
    //       + 0.32 * cos(3 * hp + pi / 30) -
    //       - 0.20 * cos(4 * hp - 63 * pi / 180)
    double hphp = hp + hp;
    double T = 1 - 0.17 * Math.cos(hp - 0.52359877559829887307710723054658)
            + 0.24 * Math.cos(hphp)
            + 0.32 * Math.cos(hphp + hp + 0.10471975511965977461542144610932)
            - 0.20 * Math.cos(hphp + hphp - 1.0995574287564276334619251841478);

    // Sh = 1 + 0.015 * Cp * T
    double Sh = 1 + 0.015 * Cp * T;

    // deltaThetaRad = (pi / 3) * e^-(36 / (5 * pi) * hp - 11)^2
    double powerBase = hp - 4.799655442984406;
    double deltaThetaRad = 1.0471975511965977461542144610932 * Math.exp(-5.25249016001879 * powerBase * powerBase);

    // Rc = 2 * sqrt((Cp^7) / (Cp^7 + 25^7))
    double Cp7 = Math.pow(Cp, 7);
    double Rc = 2 * Math.sqrt(Cp7 / (Cp7 + 6103515625.0));

    // RT = -sin(delthetarad) * Rc
    double RT = -Math.sin(deltaThetaRad) * Rc;

    // de00 = sqrt((dL / Sl)^2 + (dC / Sc)^2 + (dH / Sh)^2 + RT * (dC / Sc) * (dH / Sh))
    double dLSl = dL / Sl;
    double dCSc = dC / Sc;
    double dHSh = dH / Sh;
    return Math.sqrt(dLSl * dLSl + dCSc * dCSc + dHSh * dHSh + RT * dCSc * dHSh);
}

【问题讨论】：

你确定你需要ciede2000而不是简单的CIELAB差异吗？它的计算成本很高，并且只能为您提供“科学地”可观察到的差异。对于大多数图像处理任务而言，与 CIELAB 颜色分量的简单欧几里得距离就足够了。

标签： java algorithm optimization colors color-space

【解决方案1】：

cos 很贵，尤其是连续 4 个。您似乎正在计算 cos(na+b) 其中 b 是一个常数，n 是一个小整数。这意味着您可以预先计算 cos(b) 和 sin(b)，并在运行时只计算 cos(hp) 和 sin(hp)。你可以通过重复使用得到cos(na+b)

cos(a+b) = cos(a)*cos(b)-sin(a)*sin(b)

您将用几个sins 和coss 换取一些乘法和加法，这几乎肯定是值得的。

如果你有雄心壮志，你可以做得更好。您从atan2 间接获得hp。模式trig-function(rational-function(inverse-trig-function(x))) 经常可以替换为多项式和根的组合，它们的计算速度比三角函数更快。

我不知道pow在Java中是如何实现的，但是如果它使用日志，你最好使用Cp2=Cp*Cp;Cp4=Cp2*Cp2;Cp7=Cp4*Cp2*Cp;来获取Cp7

更新：由于我没有时间实际重写代码，因此现在变得更加投机。功率优化和三角优化其实是变相的一回事！三角优化是应用于复数的幂优化的一个版本。更何况这条线

double dH = 2 * Math.sqrt(CpProd) * Math.sin(dhp / 2);

是复数平方根运算的一部分。这使我认为，实际上可以编写大部分代码来使用复数，从而消除几乎所有的三角函数。我不知道你的复数算术如何......

【讨论】：

我已经添加了电源优化，现在功能快了1.5倍。现在我需要抓起一张纸草拟一些三角数学。 (:
余弦展开确实有点帮助，但我的复杂算术从来都不是很好，恐怕它永远也不会。 ):

【解决方案2】：

通常，任何实现此功能并存在严重速度问题的系统都不会使用随机颜色。它将做几种不同的颜色。即使是一张充满不同颜色的巨型图像，通常也只有几千种颜色。我非常强烈推荐一种缓存算法。虽然如果速度是一个问题，你应该自己总结（你只想要原语，速度）。

对于实际的颜色距离例程本身并没有太多优化，但我为这个东西编写了一个缓存系统，它的速度提高了大约 100 倍。距离惯例从压倒性的主导因素变成了昙花一现。你不应该寻求降低速度。你可能会发现一些东西。但是，请减少正确调用事物的次数。

你有两个集合输入，它产生一个集合输出，并且在很长一段时间后才会这样做。每个缓存索引 7 双。那是 14 个字节。对于 14 meg 的内存占用（左右，忽略散列或其他，可能我们说的是双倍的）。您可以存储一百万个条目，如果您有 1k 种典型的不同颜色，您将获得高达 90%s 的缓存命中率，这就足够了。如果您将初始颜色从 RGB 转换为 Lab，您甚至可以大大减少这种情况（这些转换也应该被缓存）。如果你击中 5% 的时间，你会看到加速。而且你可能会在 99% 的时间内获得成功（除非你正在做一些奇怪的事情，比如随机颜色比较）。根据我的观察，它使 CIEDE2000 花费的时间与欧几里得 RGB 几乎相同。

【讨论】：