色の距離（色差）の計算方法

はじめに

例えば以下の2つの画像，左と右の色を比べたときにどちらの方が色が近いと思いますか？

#FF0000と#400000

#FF0000と#FFBF00

単純にRGB値の距離を計算すると二つの色の距離は同じ．
しかし１つ目の画像は単に暗くなっただけなのに対して，２つ目の画像は色相が変わっているため，例えば写真のマッチングを取る場合は上の方が似た色として判別した方が都合が良くなったりする．

こうした色の距離，色差の計算について実装してみた．
ソースはJava 8を前提に書いているが，Androidでも動作する（はず）．ただし，AndroidではColorクラスが違うため適宜変更してほしい．

下準備

JavaとAndroidのColorクラスの違いや，RGB値の保存などを目的として下記のクラスを準備する．RGB値のint値を受け取るとそれぞれをdoubleの値として配列を作成し，保存する．白と黒については後で使用するためstaticな値として保持している．

DColor.java

import java.awt.*;

public class DColor {
    public static final DColor WHITE = new DColor(Color.WHITE.getRGB());
    public static final DColor BLACK = new DColor(Color.BLACK.getRGB());

    public static double getRed(int rgb){
        return ((0xff0000 & rgb) >> 16) / 255d;
    }

    public static double getGreen(int rgb){
        return ((0xff00 & rgb) >> 8) / 255d;
    }

    public static double getBlue(int rgb){
        return (0xff & rgb) / 255d;
    }

    private final double[] color;

    public DColor(int rgb){
        color = new double[]{getRed(rgb), getGreen(rgb), getBlue(rgb)};
    }

    public DColor(double r, double g, double b){
        color = new double[]{r, g, b};
    }

    public double r(){
        return color[0];
    }

    public double g(){
        return color[1];
    }

    public double b(){
        return color[2];
    }

    public double[] color(){
        return color;
    }
}

後，色差を取得するインタフェースを定義しておく．

ColorDifference.java

public interface ColorDifference {
    public double difference(DColor src, DColor dst);
}

RGB表色系でのユークリッド距離による色差の計算

まずは単純なユークリッド距離による計算クラス．色差の最大となる白と黒の距離をMAXとして，0〜1で正規化する．後述する他のクラスでも同様．

SimpleColorDifference.java

public class SimpleColorDifference implements ColorDifference {
    private static final double MAX = Math.sqrt(3);

    @Override
    public double difference(DColor src, DColor dst) {
        double rd = src.r() - dst.r();
        double gd = src.g() - dst.g();
        double bd = src.b() - dst.b();

        return Math.sqrt(rd * rd + gd * gd + bd * bd) / MAX;
    }

さて，これで最初の例を計算してみよう．

Main.java

import java.awt.*;

public class Main {
    public static void main(String[] args){
        DColor c1 = new DColor((new Color(255, 0, 0)).getRGB());
        DColor c2 = new DColor((new Color(64, 0, 0)).getRGB());
        DColor c3 = new DColor((new Color(255, 191, 0)).getRGB());

        ColorDifference cd = new SimpleColorDifference();
        System.out.println(cd.difference(c1, c2)+"\t"+cd.difference(c1, c3));
    }
}

実行結果（RGB）

0.43244667221654326 0.43244667221654326

ということで，2つの距離は等しくなった．

XYZ表色系でのユークリッド距離による色差の計算

さて，ここからはぶっちゃけ筆者がよく分かっていないのだが，色の表現としてRGB表色系ではなくXYZ表色系が用いられるらしい．他にも有名どころだとCMY表色系とかがあるが，この後にでてくる$L^* a^* b^*$表色系に繋げるためにXYZ表色系でのユークリッド距離を測ってみる．
内容は若干難読だけれどやってることは単純で単に3×3の行列をかけて変換するだけになる．

XYZDifference.java

public class XYZDifference implements ColorDifference {
    private static final double[][] m, mi;
    private static final double MAX;

    static {
        m = new double[][]{
                {0.4124f, 0.3576f, 0.1805f},
                {0.2126f, 0.7152f, 0.0722f},
                {0.0193f, 0.1192f, 0.9505f}};
        double m11=m[0][0],m12=m[0][1],m13=m[0][2];
        double m21=m[1][0],m22=m[1][1],m23=m[1][2];
        double m31=m[2][0],m32=m[2][1],m33=m[2][2];
        double d = m11*m22*m33
                + m21*m32*m13
                + m31*m12*m23
                - m11*m32*m23
                - m31*m22*m13
                - m21*m12*m33;
        mi = new double[][]{
                { (m22*m33-m23*m32)/d,(m13*m32-m12*m33)/d,(m12*m23-m13*m22)/d },
                { (m23*m31-m21*m33)/d,(m11*m33-m13*m31)/d,(m13*m21-m11*m23)/d },
                { (m21*m32-m22*m31)/d,(m12*m31-m11*m32)/d,(m11*m22-m12*m21)/d }
        };

        MAX = dif(DColor.WHITE, DColor.BLACK);
    }
    public static double getX(double r, double g, double b){
        return m[0][0] * r + m[0][1] * g + m[0][2] * b;
    }
    public static double getY(double r, double g, double b){
        return m[1][0] * r + m[1][1] * g + m[1][2] * b;
    }
    public static double getZ(double r, double g, double b){
        return m[2][0] * r + m[2][1] * g + m[2][2] * b;
    }
    public static double getR(double[] xyz){
        return mi[0][0] * xyz[0] + mi[0][1] * xyz[1] + mi[0][2] * xyz[2];
    }
    public static double getG(double[] xyz){
        return mi[1][0] * xyz[0] + mi[1][1] * xyz[1] + mi[1][2] * xyz[2];
    }
    public static double getB(double[] xyz){
        return mi[2][0] * xyz[0] + mi[2][1] * xyz[1] + mi[2][2] * xyz[2];
    }

    public static double[] rgb2xyz(DColor c){
        return new double[]{
                getX(c.r(), c.g(), c.b()),
                getY(c.r(), c.g(), c.b()),
                getZ(c.r(), c.g(), c.b())
        };
    }

    public static DColor xyz2rgb(double[] xyz){
        return new DColor(
                getR(xyz),
                getG(xyz),
                getB(xyz)
        );
    }

    @Override
    public double difference(DColor src, DColor dst) {
        return dif(src, dst) / MAX;
    }

    private static final double dif(DColor src, DColor dst){
        double[] s = rgb2xyz(src);
        double[] d = rgb2xyz(dst);
        double xd = s[0] - d[0];
        double yd = s[1] - d[1];
        double zd = s[2] - d[2];


        return Math.sqrt(xd * xd + yd * yd + zd * zd);
    }
}

rgb2xyz()を呼び出せばRGB表色系のからXYZ表色系に変換できる．使用していないがXYZ表色系からRGB表色系へ逆変換するxyz2rgb()も実装しておいた．行列の値は条件によっていろいろと変わるのでちゃんと研究する人は調べてから使って欲しい．とはいえ行列値が変わるぐらいなのでこのクラスは使い回せるはず．
色差の最大値は前述したWHITEとBLACKを突っ込んで保存している．

ちなみにJavaにはColorSpaceクラスが用意されており，

ColorSpace cs = ColorSpace.getInstance(ColorSpace.CS_CIEXYZ);
float[] xyz = cs.fromRGB(new float[]{1f, 0f, 0f});

とすればXYZ表色系に変換できるのだが，Androidでは使えないのと，実行速度が遅いので実装している．

さて，RGB表色系と同じように実行してみると下記のようになる．

実行結果（XYZ）

0.19789182003629588 0.34451910361308946

ということで2つめの色差は2つめの色差の1.8倍ぐらい，という結果が出た．

L*a*b*表色系でのユークリッド距離による色差の計算

さて，実はこういう色差については実際に知覚される色差をちゃんと求めましょう，という研究が昔からいろいろとされていて，1976年に色差を計るためのLab表色系が勧告されている．
ということでLab表色系に変換してみる．

LabDifference.java

public class LabDifference implements ColorDifference {
    private static double[] WP;
    private static final double t1 = Math.pow(6d/29, 3);
    private static final double d629 = 6d/29;
    private static final double d629_2 = d629 * d629;

    private static final double d296 = 29d/6 * 29d/6;
    private static final double d429 = 4d/29;
    private static final double d16116 = 16d/116;

    private static final double MAX;

    static {
        WP = XYZDifference.rgb2xyz(DColor.WHITE);
        MAX = dif(DColor.WHITE, DColor.BLACK);
    }

    public static double[] xyz2lab(double[] xyz){
        double[] dst = new double[3];
        double xxn = func(xyz[0] / WP[0]);
        double yyn = func(xyz[1] / WP[1]);
        double zzn = func(xyz[2] / WP[2]);
        dst[0] = 116f * yyn - 16;
        dst[1] = 500 * (xxn - yyn);
        dst[2] = 200 * (yyn - zzn);

        return dst;
    }

    public static double[] lab2xyz(double[] lab){
        double l = lab[0];
        double a = lab[1];
        double b = lab[2] ;
        double fy = (l+16)/116;
        double fx = fy + a / 500;
        double fz = fy - b / 200;

        return new double[]{
                fx > d629 ? WP[0] * fx * fx * fx : (fx - d16116) * 3 * d629_2 * WP[0],
                fy > d629 ? WP[1] * fy * fy * fy : (fy - d16116) * 3 * d629_2 * WP[1],
                fx > d629 ? WP[2] * fz * fz * fz : (fz - d16116) * 3 * d629_2 * WP[2]
        };
    }

    public static double func(double t){
        if(t > t1){
            return Math.pow(t, 1d/3);
        } else {
            return 1d/3 * d296 * t + d429;
        }
    }

    public static double[] rgb2lab(DColor c){
        return xyz2lab(XYZDifference.rgb2xyz(c));
    }

    public static DColor lab2rgb(double[] lab){
        return XYZDifference.xyz2rgb(lab2xyz(lab));
    }

    @Override
    public double difference(DColor src, DColor dst) {
        return dif(src, dst) / MAX;
    }

    private static final double dif(DColor src, DColor dst){
        double[] s = rgb2lab(src);
        double[] d = rgb2lab(dst);
        double xd = s[0] - d[0];
        double yd = s[1] - d[1];
        double zd = s[2] - d[2];

        return Math.sqrt(xd * xd + yd * yd + zd * zd);
    }
}

ここまで来ると門外漢にはさっぱり分からないので詳細は専門家に聞いて欲しい．ググって出てきた数式を単純に実装している．一応，実装上の工夫として固定値をあらかじめ計算しているが，Javaはコンパイル時に最適化されるはずなので速度的には意味は無く，単に実装上の誤りを減らすために行っている．

さて，Lab表色系で色差を計算すると

実行結果（Lab)

0.47194631171156237 0.9645251064099958

という結果が出た．大体2倍ぐらいの違い，というわけだ．

CIEDE2000による色差の計算

さて，本番．KONICA MINOLTAのこのページによると，Labによる表色系でもまだ問題があり，CIEDE2000というのが出来たのだが，CIEDE2000では

CIE Labの色空間上での人の目の色識別域の特徴，１）彩度依存性，２）色相依存性，および，３）明度依存性を考慮した計算指揮となっています．

だそうだ．

Lab立体色空間上で真円もしくは方形状であったのに対し、CIE2000色差式での色差ΔE00は、彩度方向を長軸とする楕円形状となり、人の目の色識別域の形状に近似しているのが特長

なんだそうで，それを地道に計算することになる．「彩度が低いと真円に近い楕円」となり「彩度が高いと彩度方向にのびた楕円」となる．つまり，人間は彩度が高いほど色差を感じにくくなる，ということらしい．確かにそんな気はする．ということで，実装してみた．
数式はWikipediaを参考にしている．他の文献でもそうなのだが，Degreeとして角度を扱っているため，全てRadiusに変換する辺りで注意が必要．（文献によっては間違えてたりもする）

CIE2kDifference.java

public class CIE2kDifference implements ColorDifference {

    private static final double v25_7 = Math.pow(25, 7);
    private static final double d6 = Math.toRadians(6);
    private static final double d25 = Math.toRadians(25);
    private static final double d30 = Math.toRadians(30);
    private static final double d60 = Math.toRadians(60);
    private static final double d63 = Math.toRadians(63);
    private static final double d275 = Math.toRadians(275);
    private static final double kl = 1;
    private static final double kc = 1;
    private static final double kh = 1;
    private static final double MAX;

    static {
        MAX = dif(DColor.WHITE, DColor.BLACK);
    }

    @Override
    public double difference(DColor src, DColor dst) {
        return dif(src, dst) / MAX;
    }

    public static double dif(DColor src, DColor dst) {
        double[] src_lab = LabDifference.rgb2lab(src);
        double l1 = src_lab[0];
        double a1 = src_lab[1];
        double b1 = src_lab[2];

        double[] dst_lab = LabDifference.rgb2lab(dst);
        double l2 = dst_lab[0];
        double a2 = dst_lab[1];
        double b2 = dst_lab[2];

        return dif(l1, a1, b1, l2, a2, b2);
    }

    public static double dif(double l1, double a1, double b1, double l2, double a2, double b2){
        double dld = l2 - l1;
        double lb = (l1 + l2) / 2;

        double cs1 = Math.hypot(a1, b1);
        double cs2 = Math.hypot(a2, b2);
        double cb = (cs1 + cs2) / 2;
        double cb7 = Math.pow(cb, 7);
        double ad1 = a1 + a1 / 2 * (1 - Math.sqrt(cb7 / (cb7 + v25_7)));
        double ad2 = a2 + a2 / 2 * (1 - Math.sqrt(cb7 / (cb7 + v25_7)));

        double cd1 = Math.hypot(ad1, b1);
        double cd2 = Math.hypot(ad2, b2);
        double cbd = (cd1 + cd2) / 2;
        double cbd7 = Math.pow(cbd, 7);

        double dcd = (cd2 - cd1);
        double hd1 = b1 == 0 && ad1 == 0 ? 0 : Math.atan2(b1, ad1);
        if(hd1 < 0){
            hd1 += Math.PI * 2;
        }
        double hd2 = b2 == 0 && ad2 == 0 ? 0 : Math.atan2(b2, ad2);
        if(hd2 < 0){
            hd2 += Math.PI * 2;
        }

        double dhd = hd2 - hd1;
        if(cd1 * cd2 == 0){
            dhd = 0;
        } else if(Math.abs(hd1 - hd2) > Math.PI) {
            if(hd2 <= hd1){
                dhd += Math.PI * 2;
            } else {
                dhd -= Math.PI * 2;
            }
        }


        double dhhd = 2 * Math.sqrt(cd1 * cd2) * Math.sin(dhd / 2);
        double hhbd = 0;
        if(cd1 * cd2 != 0){
            hhbd = Math.abs(hd1 - hd2) > Math.PI ? ( hd1 + hd2 + Math.PI * 2) / 2 : (hd1 + hd2) / 2;
        }

        double tt = 1
                - 0.17 * Math.cos(hhbd - d30)
                + 0.24 * Math.cos(2 * hhbd)
                + 0.32 * Math.cos(3 * hhbd + d6)
                - 0.20 * Math.cos(4 * hhbd - d63);
        double lb50_2 = Math.pow(lb - 50, 2);
        double ssl = 1 + (0.015 * lb50_2) / Math.sqrt(20 + lb50_2);
        double ssc = 1 + 0.045 * cbd;
        double ssh = 1 + 0.015 * cbd * tt;
        double rrt = -2d * Math.sqrt(cbd7 / (cbd7 +v25_7)) * Math.sin(d60 * Math.exp(- Math.pow((hhbd - d275)/ d25, 2)));
        double de = Math.pow(dld / (kl * ssl), 2)
                + Math.pow(dcd / (kc * ssc), 2)
                + Math.pow(dhhd / (kh * ssh), 2)
                + rrt * (dcd / (kc * ssc)) * (dhhd / (kh * ssh));

        return Math.sqrt(de);
    }
}

さて，もはやさっぱり分からないけれど，一応これで合っているらしい．この辺りの文献と比較して出力される値が同じであることは確認している．

実行結果は次の通り．

実行結果（CIE2000）

0.24093258957174238 0.5480665244873609

ということで，2.2倍ぐらいらしい．（もはや値はどうでも良くなってきたが）

おまけ（CIEDE2000の高速化）

CIEDE2000は見た目の通り凄く重い処理なので，できるだけ軽くする．とはいっても最近のJavaは最適化されるのでそもそも速く，ほとんど出来るところが無かった．
Androidだとまだまだできるところはありそうなので後日やってみる．

FastCIE2kDif.java

public class FastCIE2kDif implements ColorDifference {

    private static final double pi_2 = Math.PI * 2;
    private static final double v25_7 = Math.pow(25, 7);
    private static final double d6 = Math.toRadians(6);
    private static final double d25 = Math.toRadians(25);
    private static final double d30 = Math.toRadians(30);
    private static final double d60 = Math.toRadians(60);
    private static final double d63 = Math.toRadians(63);
    private static final double d275 = Math.toRadians(275);
    private static final double kl = 1;
    private static final double kc = 1;
    private static final double kh = 1;
    private static final double MAX;

    static {
        MAX = dif(DColor.WHITE, DColor.BLACK);
    }

    public static final double hypot(double a, double b){
        return Math.sqrt(a * a + b * b);
    }

    public static final double pow2(double a){
        return a * a;
    }

    public static final double pow7(double a){
        return a * a * a * a * a * a * a;
    }

    @Override
    public double difference(DColor src, DColor dst) {
        return dif(src, dst) / MAX;
    }

    public static double dif(DColor src, DColor dst) {
        double[] src_lab = LabDifference.rgb2lab(src);
        double l1 = src_lab[0];
        double a1 = src_lab[1];
        double b1 = src_lab[2];

        double[] dst_lab = LabDifference.rgb2lab(dst);
        double l2 = dst_lab[0];
        double a2 = dst_lab[1];
        double b2 = dst_lab[2];

        return dif(l1, a1, b1, l2, a2, b2);
    }

    public static double dif(double l1, double a1, double b1, double l2, double a2, double b2){
        double dld = l2 - l1;
        double lb = (l1 + l2) / 2;

        double cs1 = hypot(a1, b1);
        double cs2 = hypot(a2, b2);
        double cb = (cs1 + cs2) / 2;
        double cb7 = pow7(cb);
        double ad1 = a1 + a1 / 2 * (1 - Math.sqrt(cb7 / (cb7 + v25_7)));
        double ad2 = a2 + a2 / 2 * (1 - Math.sqrt(cb7 / (cb7 + v25_7)));

        double cd1 = hypot(ad1, b1);
        double cd2 = hypot(ad2, b2);
        double cbd = (cd1 + cd2) / 2;
        double cbd7 = pow7(cbd);

        double dcd = (cd2 - cd1);
        double hd1 = b1 == 0 && ad1 == 0 ? 0 : Math.atan2(b1, ad1);
        if(hd1 < 0){
            hd1 += pi_2;
        }
        double hd2 = b2 == 0 && ad2 == 0 ? 0 : Math.atan2(b2, ad2);
        if(hd2 < 0){
            hd2 += pi_2;
        }

        double dhd = hd2 - hd1;
        if(cd1 * cd2 == 0){
            dhd = 0;
        } else if(Math.abs(hd1 - hd2) > Math.PI) {
            if(hd2 <= hd1){
                dhd += pi_2;
            } else {
                dhd -= pi_2;
            }
        }


        double dhhd = 2 * Math.sqrt(cd1 * cd2) * Math.sin(dhd / 2);
        double hhbd = 0;
        if(cd1 * cd2 != 0){
            hhbd = Math.abs(hd1 - hd2) > Math.PI ? ( hd1 + hd2 + pi_2) / 2 : (hd1 + hd2) / 2;
        }

        double tt = 1
                - 0.17 * Math.cos(hhbd - d30)
                + 0.24 * Math.cos(2 * hhbd)
                + 0.32 * Math.cos(3 * hhbd + d6)
                - 0.20 * Math.cos(4 * hhbd - d63);
        double lb50_2 = pow2(lb - 50);
        double ssl = 1 + (0.015 * lb50_2) / Math.sqrt(20 + lb50_2);
        double ssc = 1 + 0.045 * cbd;
        double ssh = 1 + 0.015 * cbd * tt;
        double rrt = -2d * Math.sqrt(cbd7 / (cbd7 +v25_7)) * Math.sin(d60 * Math.exp(- pow2((hhbd - d275)/ d25)));
        double de = pow2(dld / (kl * ssl))
                + pow2(dcd / (kc * ssc))
                + pow2(dhhd / (kh * ssh))
                + rrt * (dcd / (kc * ssc)) * (dhhd / (kh * ssh));

        return Math.sqrt(de);
    }
}

やったことは下記の通り．

• Math.hypot()が凄く遅いのでベタで書き直した．バッファオーバーフローした場合にちゃんと計算できないが今回の例では大きな値を使わないので問題無い．
• Math.pow()では2乗と7乗しか使ってないのでそれぞれメソッドを用意して掛け算を行った．ただ，2乗に関しては効果無し，7乗の内部を2->4->8->7乗の順に計算させたりいろいろやったが単純に7回掛けるのと変化無し．

他にもcosやatanを近似計算させたりしてみたが全然効果が出なかった．角度計算をint化して最後にRadius変換する方法等も考えたが効果が出そうに無いのでやっていない．

やることリスト

GitHubに上げたりAndroid用ライブラリを作る予定．