RGB→LAB変換に見るPythonとC++の速度差

PythonからC++のクラスを使う方法は前回学んだが、
dllのAPIを大量に叩いた時に実用に耐えるのかは気になります。
http://qiita.com/tibigame/items/499a81f0e1cb47962b90
今回はRGB→LABの変換とRGB2色からCIEDE2000色差を求める処理で速度を計る。

なお色変換のJavascriptによる実装はこちらを参照していただきたい。
http://qiita.com/tibigame/items/40ab217c863a20cdb264
ちなみにES6記法でエレガントに書いているためか相当速くて優秀だった。

Pythonによる色変換の実装

実験なのでコア部分だけを示す。
俺はJavascript配列のmapメソッドを使いたいんだ的な実装である。
おかげでコードの構造はほぼ同じになった。

Color.py

import math
import numpy as np
import enum

ColorType = enum.Enum("ColorType", "RGB HEX LAB LCH")
class Color:
    def __lab_to_rgb(self, lab):
        l, a, b = lab[0], lab[1], lab[2]
        y = (l * 100) / 903.3 if l <= 8 else 100 * ((l + 16) / 116) ** 3
        y2 = (7.787 * (y / 100)) + (16 / 116) if l <= 8 else np.cbrt(y / 100.0)
        x = (95.047 * ((a / 500) + y2 - (16 / 116))) / 7.787 if a / 500.0 + y2 <= 0.2069 else 95.047 * ((a / 500) + y2) ** 3
        z = (108.883 * (y2 - (b / 200) - (16 / 116))) / 7.787 if y2 - b / 200.0 <= 0.2059 else 108.883 * (y2 - (b / 200)) ** 3
        ndXYZ = np.matmul(np.array(
        [[0.032406, -0.015372, -0.004986],
        [-0.009689, 0.018758, 0.000415],
        [0.000557, -0.002040, 0.010570]]), np.array([x, y, z]).T)
        ndCond = ndXYZ > 0.0031308
        ndTrue = (1.055 * np.power(ndXYZ, 1.0 / 2.4)) - 0.055
        ndFalse = ndXYZ * 12.92
        return list(np.minimum(np.maximum(0, np.where(ndCond, ndTrue, ndFalse)), 1.0) * 255.0)

colorDifferenceType = enum.Enum("colorDifferenceType", "EuclidLAB EuclidRGB CIEDE2000")
def colorDifference(Color1: Color, Color2: Color, colorDifferenceType = colorDifferenceType.EuclidLAB):
    if colorDifferenceType == colorDifferenceType.EuclidLAB:
        Color1.Calc(ColorType.LAB)
        Color2.Calc(ColorType.LAB)
        return np.linalg.norm(np.array(Color1.LAB) - np.array(Color2.LAB))
    elif colorDifferenceType == colorDifferenceType.EuclidRGB:
        Color1.Calc(ColorType.RGB)
        Color2.Calc(ColorType.RGB)
        return np.linalg.norm(np.array(Color1.RGB) - np.array(Color2.RGB))
    elif colorDifferenceType == colorDifferenceType.CIEDE2000:
        tolerance_zero = (lambda x: np.abs(x) < 0.00000001)
        cosd = (lambda deg: np.cos(np.deg2rad(deg)))
        sind = (lambda deg: np.sin(np.deg2rad(deg)))
        reg_fqatan = (lambda deg: deg if deg >= 0 else deg + 360.0)
        fqatan = (lambda x, y: reg_fqatan(np.rad2deg(np.arctan2(x, y))))
        f7 = (lambda x: np.power(x / 25.0, 3.5) if x < 1.0 else 1.0 / np.sqrt(1.0 + np.power(25.0 / x, 7.0)))

        Color1.Calc(ColorType.LAB)
        Color2.Calc(ColorType.LAB)
        L1, a1, b1 = Color1.LAB
        L2, a2, b2 = Color2.LAB
        epsilon = 0.000000001
        C1ab = np.linalg.norm(np.array([a1, b1]))
        C2ab = np.linalg.norm(np.array([a2, b2]))
        G = 0.5 * (1.0 - f7((C1ab+C2ab) / 2.0)) + 1.0
        a1_ = G * a1
        a2_ = G * a2
        C1_ = np.linalg.norm(np.array([a1_, b1]))
        C2_ = np.linalg.norm(np.array([a2_, b2]))
        h1_ = 0.0 if tolerance_zero(a1_) and tolerance_zero(b1) else fqatan(b1, a1_)
        h2_ = 0.0 if tolerance_zero(a2_) and tolerance_zero(b2) else fqatan(b2, a2_)
        dL_ = L2 - L1
        dC_ = C2_ - C1_
        C12 = C1_ * C2_
        if tolerance_zero(C12):
            dh_ = 0.0
        else:
            tmp = h2_ - h1_
            if np.abs(tmp) <= 180.0 + epsilon:
                dh_ = tmp
            elif tmp > 180.0:
                dh_ = tmp - 360.0
            elif tmp < -180.0:
                dh_ = tmp + 360.0
        dH_ = 2.0 * np.sqrt(C12) * sind(dh_ / 2.0)

        if tolerance_zero(C12):
            h_ = h1_ + h2_
        else:
            tmp2 = h1_ + h2_
            if np.abs(h1_ - h2_) <= 180.0 + epsilon:
                h_ = tmp2 / 2.0
            elif tmp2 < 360.0:
                h_ = (tmp2 + 360.0) / 2.0
            elif tmp2 >= 360.0:
                h_ = (tmp2 - 360.0) / 2.0
        L_ = (L1 + L2) / 2.0
        C_ = (C1_ + C2_) / 2.0
        L_2 = (L_ - 50.0) * (L_ - 50.0)
        SL = 1.0 + 0.015 * L_2 / np.sqrt(20.0 + L_2)
        SC = 1.0 + 0.045 * C_
        SH = 1.0 + 0.015 * C_ * (1.0 - 0.17 * cosd(h_ - 30.0) + 0.24 * cosd(2.0 * h_) + 0.32 * cosd(3.0 * h_ + 6.0) - 0.2 * cosd(4.0 * h_ -63.0))
        LP = dL_ / SL
        CP = dC_ / SC
        HP = dH_ / SH
        expInner = (h_ - 275.0) / 25.0
        RT = - sind(2.0 * 30.0 * np.exp(- expInner * expInner)) * 2.0 * f7(C_)
        return np.sqrt(LP * LP + CP * CP + HP * HP + RT * CP * HP)

C++による実装

こちらもラムダ式を使って脳筋でJavascriptのコードをそのまま実装した。

Colordll.cpp

class Color {
public:
    std::valarray<double> rgb_to_lab() {
        std::valarray<double> rgbValarray = { RGB[0] / 255.0, RGB[1] / 255.0, RGB[2] / 255.0 };
        std::valarray<double> rgbD = rgbValarray.apply([](double d) { return d > 0.04045 ? std::pow((d + 0.055) / 1.055, 2.4) : d / 12.92; });
        std::valarray<double> xyz = { 
            ((rgbD[0] * 0.4124) + (rgbD[1] * 0.3576) + (rgbD[2] * 0.1805)) / 0.95047,
            ((rgbD[0] * 0.2126) + (rgbD[1] * 0.7152) + (rgbD[2] * 0.0722)),
            ((rgbD[0] * 0.0193) + (rgbD[1] * 0.1192) + (rgbD[2] * 0.9505)) / 1.08883 };
        std::valarray<double> xyzD = xyz.apply([](double d) { return d > 0.008856 ? std::cbrt(d) : 7.787 * d + 16 / 116; });
        std::valarray<double> result = { (116.0 * xyzD[1]) - 16.0, 500.0 * (xyzD[0] - xyzD[1]), 200.0 * (xyzD[1] - xyzD[2]) };
        return result;
    }
}

double colorDiffCIEDE2000Class(Color A, Color B) {
    auto cosd = [](double deg) { return std::cos(deg * M_PI / 180); };
    auto sind = [](double deg) { return std::sin(deg * M_PI / 180); };
    auto reg_fqatan = [](double deg) { return deg >= 0 ? deg : deg + 360.0; };
    auto fqatan = [reg_fqatan](double x, double y) { return reg_fqatan(std::atan2(x, y) * 180.0 / M_PI); };
    auto f7 = [](double x) { return x < 1.0 ? std::pow(x / 25.0, 3.5) : 1.0 / std::sqrt(1.0 + std::pow(25.0 / x, 7)); };

    const std::valarray<double> LAB1 = A.getDouble(SPACE_LAB);
    const std::valarray<double> LAB2 = B.getDouble(SPACE_LAB);
    const double L1 = LAB1[0];
    const double a1 = LAB1[1];
    const double b1 = LAB1[2];
    const double L2 = LAB2[0];
    const double a2 = LAB2[1];
    const double b2 = LAB2[2];

    const double C1ab = std::hypot(a1, b1);
    const double C2ab = std::hypot(a2, b2);
    const double Cab = (C1ab + C2ab) / 2.0;
    const double G = 0.5 * (1.0 - f7(Cab));
    const double a1_ = (1.0 + G) * a1;
    const double a2_ = (1.0 + G) * a2;
    const double C1_ = std::hypot(a1_, b1);
    const double C2_ = std::hypot(a2_, b2);
    const double h1_ = fqatan(b1, a1_);
    const double h2_ = fqatan(b2, a2_);
    const double dL_ = L2 - L1;
    const double dC_ = C2_ - C1_;
    const double C12 = C1_ * C2_;

    double dh_;
    const double tmp = h2_ - h1_;
    if (std::abs(tmp) <= 180.0) {
        dh_ = tmp;
    }
    else if (tmp > 180.0) {
        dh_ = tmp - 360.0;
    }
    else {
        dh_ = tmp + 360.0;
    }

    const double dH_ = 2.0 * std::sqrt(C12) * sind(dh_ / 2.0);
    const double L_ = (L1 + L2) / 2.0;
    const double C_ = (C1_ + C2_) / 2.0;

    double h_;
    const double tmp1 = std::abs(h1_ - h2_);
    const double tmp2 = h1_ + h2_;
    if (tmp1 <= 180.0) {
        h_ = tmp2 / 2.0;
    }
    else if (tmp2 < 360.0) {
        h_ = (tmp2 + 360.0) / 2.0;
    }
    else {
        h_ = (tmp2 - 360.0) / 2.0;
    }

    const double T = 1.0 - 0.17 * cosd(h_ - 30.0) + 0.24 * cosd(2.0 * h_)
        + 0.32 * cosd(3.0 * h_ + 6.0) - 0.2 * cosd(4.0 * h_ - 63.0);
    const double dTh = 30.0 * std::exp(-std::pow((h_ - 275.0) / 25.0, 2.0));
    const double L_2 = (L_ - 50.0) * (L_ - 50.0);
    const double RC = 2.0 * f7(C_);
    const double SL = 1.0 + 0.015 * L_2 / std::sqrt(20.0 + L_2);
    const double SC = 1.0 + 0.045 * C_;
    const double SH = 1.0 + 0.015 * C_ * T;
    const double RT = -sind(2.0 * dTh) * RC;
    const double kL = 1.0;
    const double kC = 1.0;
    const double kH = 1.0;
    const double LP = dL_ / (kL * SL);
    const double CP = dC_ / (kC * SC);
    const double HP = dH_ / (kH * SH);
    return std::sqrt(LP * LP + CP * CP + HP * HP + RT * CP * HP);
}

Pythonでの計測用コード

テストはループを256*256*12=786432回まわすことで行う。(6コア12スレッドマシンでの計測の都合)
[loopFunc(i, j, k) for k in range(12) for i in range(256) for j in range(256)]
のようなリスト内包形式は優位な時間差がなかったので可読性からforループを採用した。

test.py

import ColorLib
from Color import Color, ColorType, colorDifference, colorDifferenceType, rgb_to_lab
import time

# 時間計測のベース関数
def test_base(func, text):
    starttime = time.time()
    func()
    endtime = time.time()
    interval = starttime - endtime
    print (str(-interval) + "秒 " + text)
# テストされる関数
def test_python_normal():
    for i in range(256):
        for j in range(256):
            for k in range(12):
                c = Color([[k, j, i], 1.0], ColorType.RGB)
                c.Calc(ColorType.LAB)
def test_Cdll():
    for i in range(256):
        for j in range(256):
            for k in range(12):
                y = ColorLib.Color([k, j, i], ColorLib.SPACE_RGB, 1.0)
                y.rgb_to_lab()
def test_Python_CIEDE2000():
    for i in range(256):
        for j in range(256):
            for k in range(12):
                c = Color([[k, j, i], 1.0], ColorType.RGB)
                d = Color([[j, i, k], 1.0], ColorType.RGB)
                c.Calc(ColorType.LAB)
                d.Calc(ColorType.LAB)
                colorDifference(c, d, colorDifferenceType.CIEDE2000)
# このコードはクラス生成と計算を2回ずつ、計算を1回で合計5回のAPIをループごとに呼んでいる
def test_Cdll_CIEDE2000():
    for i in range(256):
        for j in range(256):
            for k in range(12):
                y = ColorLib.Color([k, j, i], ColorLib.SPACE_RGB, 1.0)
                C1 = y.rgb_to_lab()
                x = ColorLib.Color([j, i, k], ColorLib.SPACE_RGB, 1.0)
                C2 = x.rgb_to_lab()
                ColorLib.colorDiffCIEDE2000([C1[0], C1[1], C1[2], C2[0], C2[1], C2[2]])
# 仮に1回だとどうなるかも試してみる
def test_Cdll_CIEDE2000_Func():
    for i in range(256):
        for j in range(256):
            for k in range(12):
                ColorLib.colorDiffCIEDE2000Func([k,i,j,j,i,k])
# APIを1回しか呼ばない=全部C++でやる
def test_rgb_to_lab_pureC():
    c = ColorLib.rgb_to_lab_pureC()

# テストを実行する
test_base(test_python_normal, "(Python)")
test_base(test_Cdll, "(Python C++dll)")
test_base(test_Python_CIEDE2000, "(Python CIEDE2000)")
test_base(test_Cdll_CIEDE2000, "(Python C++dll CIEDE2000)")
test_base(test_Cdll_CIEDE2000_Func, "(Python C++dll CIEDE2000_Func)")
test_base(test_rgb_to_lab_pureC, "(C++)")

Pythonでのなんちゃってマルチスレッド

Pythonは簡単にマルチプロセスで計算できると聞いてやってみた。
Poolしてmapするだけとか超簡単だ。

testMulti.py

import ColorLib
from Color import Color, ColorType, colorDifference, colorDifferenceType
import time

from multiprocessing import Pool
from multiprocessing import Process, freeze_support

def test_python_normal(k):
    for i in range(256):
        for j in range(256):
            c = Color([[k, j, i], 1.0], ColorType.RGB)
            c.Calc(ColorType.LAB)
def test_python_normal_dll(k):
    for i in range(256):
        for j in range(256):
            y = ColorLib.Color([k, j, i], ColorLib.SPACE_RGB, 0.4281142)
            y.rgb_to_lab()

if __name__ == '__main__':
    freeze_support()
    starttime = time.time()
    p = Pool()
    p.map(test_python_normal, range(12))
    endtime = time.time()
    interval = starttime - endtime
    print (str(-interval) + "秒 (Python) マルチスレッド")

    starttime = time.time()
    p = Pool()
    p.map(test_python_normal_dll, range(12))
    endtime = time.time()
    interval = starttime - endtime
    print (str(-interval) + "秒 (C++dll) マルチスレッド")

結果

実装方式	RGB→Lab関数	RGB→Labクラス	CIEDE2000関数	CIEDE2000クラス
C++	0.82	0.86		2.86
Javascript	2.57		4.15
Python(dll)MP		1.15	1.23	5.39
Python MP		4.7		17.5
Python(dll)		2.56	5.39	10.41
Python	18.36	25.31		117.28

C++はクラス構造とか十分に最適化の余地はある。
なによりシングルスレッドだし。
Javascriptは型変換の入らない計算ばかりのためか大健闘している。
Pythonでのなんちゃってマルチプロセスは十分に機能している。
素のPythonではコア数倍を超えてスレッド数倍近くまで成績を伸ばしていて費用対効果抜群だ。
C++dllコールはせいぜい倍程度でマルチプロセスの効果は少ない。
これはdllコールのコストが高いからと推測され、
実際クラスと関数の差を見るとdllコールを1/5に減らした関数版だと
シングル・マルチともに大幅に改善されていることがわかる。

結局PythonからC++のdllを叩くのが勝手が良さそうだ。
マルチプロセスでブーストすれば並のC++プログラム相手でも戦えるし、
いざとなったらdllのコールを減らすための関数をC++側に実装すればいい。

ちなみにCythonでCIEDE2000関数を実装したところ25秒ほどかかった。
dll版の5倍だ。
これは純粋にコンパイラの差に思える。
AVX2オプションなどでの最適化は想像以上に大きい。
pow関数をstd::powからAVX2でSIMD実装されたものに置き換えるだけで10%ぐらい速くなったりもするし。