はじめに

Kmeansの拡張でkを指定しないXmeansというのがありますが
本記事では画像の減色する際に類似度を指定して、色数を決めてみたいと思います。
普通にXmeansでやろうもんなら遅くて使いもんにならないだろうし、クラスタ数ではなく類似度で指定したいなと思ったので、考えた次第です。

アルゴリズム概要

画像を縮小します。
kMin から kMax まで1ずつ増やして Kmeans 法を行います。
Kmeans の前後の画像で類似度を計算し、それが類似度が指定されたものより大きければそれを採用します。

類似度について

$d({\bf x}, {\bf y}) / \sqrt{256^3\cdot 3}<\epsilon$を満たすとき、画素の色が一致するとします。
$d:{\mathbb R}^2\times{\mathbb R}^2\rightarrow{\mathbb R}$ はユークリッド距離関数とします。
画像の類似度は、一致する画素の数を全ての画素の数で割ったものとします。

実験

実際に動かしたものを載せます。
変数は

const [kMin, kMax, colorEpsilon, size, bitDepth, attempts] = [2, 16, 0.08, 32, 4, 2];

のように固定します。similarity(類似度)を変更したものが以下です。
処理時間はほぼ変わらないので記載していません。全体で100msで複数回のKmeansで20msぐらいです。

減色後の画像	類似度	色数
	0.8	4
	0.85	5
	0.9	5
	0.95	7
	0.99	12

類似度が大きくすれば、おおむね色数が増えていることが確認できました。

ソース

    /**
     * Determine the number of colors based on similarity
     * @param {Mat} src source mat 
     * @param {Mat} dst destination mat
     * @param {number} kMin cluster min count
     * @param {number} kMax cluster max count
     * @param {number} similarity similarity (ex. 0.95)
     * @param {number} colorEpsilon color tolerance (ex. 0.08)  
     * @param {Size} dsize size in resizing
     * @param {number} interpolation interpolation
     * @param {number} bitDepth bit depth [1, 8]
     * @param {number} attempts number of attempts 
     * @param {number} maxCount the maximum number of iterations/elements
     * @param {number} epsilon the desired accuracy
     * @returns {Mat} destination mat
     */
    static fastXmeans(src, dst, kMin, kMax, similarity, colorEpsilon, 
        dsize, interpolation, bitDepth, attempts, maxCount = 1e4, epsilon = 1e-4) {
        // Resize
        const resized = new cv.Mat();
        cv.resize(src, resized, dsize, 0, 0, interpolation);

        // Initialize
        dst = dst || new cv.Mat(src.size(), src.type());
        const criteria = new cv.TermCriteria(cv.TermCriteria_EPS + cv.TermCriteria_MAX_ITER, maxCount, epsilon);
        const flags = cv.KMEANS_PP_CENTERS;

        // Create samples
        const samples = new cv.Mat(resized.rows * resized.cols, 3, cv.CV_32F);
        for(let y = 0; y < resized.rows; y += 1) {
            for(let x = 0; x < resized.cols; x += 1) {
                const idx = y * resized.cols + x;
                for(let z = 0; z < 3; z += 1) {                    
                    samples.floatPtr(idx)[z] = resized.data[idx * 4 + z];
                }
            }
        }

        // Decide k(cluster count)
        let startDate = new Date();
        let x = -1;
        const labelsArray = [], ui8CentersArray = [];
        for(let k = kMin; k <= kMax; k += 1) {
            // K-Means method            
            const [labels, centers] = [new cv.Mat(), new cv.Mat()];
            cv.kmeans(samples, k, labels, criteria, attempts, flags, centers);
            labelsArray.push(labels);

            // Convert Float to UInt8
            const floatToUint8 = x => Math.min(Math.max(Math.round(x), 0), 255);
            const ui8Centers = new Uint8Array(k * 3);
            for(let c = 0; c < k; c += 1) {
                for(let z = 0; z < 3; z += 1) {
                    ui8Centers[c * 3 + z] = floatToUint8(centers.floatAt(c, z));
                }
            }
            centers.delete();
            ui8CentersArray.push(ui8Centers);
            
            // Create kdst from labels and centers
            const kdst = new cv.Mat(resized.size(), resized.type());
            for(let y = 0; y < kdst.rows; y += 1) {
                for(let x = 0; x < kdst.cols; x += 1) {
                    const p = (y * kdst.cols + x) * 4;
                    const c = labels.data[p];
                    for(let z = 0; z < 3; z += 1) {                    
                        kdst.data[p + z] = ui8Centers[c * 3 + z];
                    }
                    kdst.data[p + 3] = resized.data[p + 3];
                }
            }
            // Compute sim            
            const sim = computeSimilarity(resized, kdst, colorEpsilon);
            if(sim > similarity) {
                x = k;
                break;
            }
        }
        if(x < 0) { x = kMax; } // not found.
        const k = x;
        const labels = labelsArray[x - kMin];
        const ui8Centers = ui8CentersArray[x - kMin];
        console.log('kmeans:', new Date() - startDate, 'ms');
        
        // Create color map
        const mem = 2 ** bitDepth;
        const quant = (mem - 1) / (256 - 1);
        const digit = Math.ceil(Math.log10(mem - 1));
        const colorMap = {};
        for(let y = 0; y < resized.rows; y += 1) {
            for(let x = 0; x < resized.cols; x += 1) {
                const p = (y * resized.cols + x) * 4;
                const key = createColorKey(resized.data, p);
                if(typeof colorMap[key] === 'undefined') {
                    colorMap[key] = labels.data[p];
                }
            }
        }        
        // Release
        resized.delete();
        labelsArray.forEach(e => { e.delete(); });

        // Create dst from ui8Centers
        startDate = new Date();
        for(let y = 0; y < dst.rows; y += 1) {
            for(let x = 0; x < dst.cols; x += 1) {
                const p = (y * dst.cols + x) * 4;
                const key = createColorKey(src.data, p);
                // Find nearest color
                if(typeof colorMap[key] === 'undefined') {
                    let [minDist2, minC] = [Number.MAX_VALUE, -1];
                    for(let c = 0; c < k; c += 1) {
                        // Compute distance squared in color space
                        let dist2 = 0;
                        for(let z = 0; z < 3; z += 1) {
                            dist2 += (src.data[p + z] - ui8Centers[c * 3 + z]) ** 2;
                        }
                        if(dist2 < minDist2) {
                            minDist2 = dist2;
                            minC = c;
                        }
                    }
                    colorMap[key] = minC;
                }          
                // Set color
                const c = colorMap[key];
                for(let z = 0; z < 3; z += 1) {
                    dst.data[p + z] = ui8Centers[c * 3 + z];
                }
                dst.data[p + 3] = src.data[p + 3];
            }
        }
        console.log('after kmeans:', new Date() - startDate, 'ms');
        return [ dst, k ];

        /**
         * Create color key
         * @param {Uint8Array} data data 
         * @param {number} p pixel index 
         * @returns {string} color key
         */
        function createColorKey(data, p) {
            let key = '';
            for(let z = 0; z < 3; z += 1) {
                const c = Math.round(data[p + z] * quant); // Quantization
                key += c.toString().padStart(digit, '0');
            }
            return key;     
        }
        /**
         * Compute image similarity
         * @param {Mat} src source mat
         * @param {Mat} dst destination mat 
         * @param {number} epsilon similarity epsilon ex. 0.08
         * @returns {number} image similarity
         */
        function computeSimilarity(src, dst, epsilon) {
            const sqrt3 = 3 ** 0.5;
            let diffCount = 0;
            for(let p = 0; p < src.data.length; p += 4) {
                const srcColor = getColor(src.data, p);
                const dstColor = getColor(dst.data, p);
                const diff = diffColor(srcColor, dstColor);
                if(diff / (255 * sqrt3) > epsilon) {
                    diffCount++;
                } 
            }
            return 1 - diffCount / (src.cols * src.rows);

            function diffColor(src, dst) {
                return ((src[0] - dst[0]) ** 2 + (src[1] - dst[1]) ** 2 + (src[2] - dst[2]) ** 2) ** 0.5;
            }
            function getColor(data, p) {
                return [ data[p], data[p + 1], data[p + 2] ];
            }  
        }              
    }

関数の呼び出し方

// Create src and dst
const src = cv.imread('src-canvas');
const fx = fastXmeans;
const [kMin, kMax, similarity, colorEpsilon, size, bitDepth, attempts] = [2, 16, 0.9, 0.08, 32, 4, 2];
let [dst, k] = OpenCVWrapper.fastXmeans(src, null, kMin, kMax, similarity, colorEpsilon, 
                new cv.Size(size, size), cv.INTER_LINEAR, bitDepth, attempts);
            
// Show dst
cv.imshow('dst-canvas', dst);

// Remove src and dst
src.delete();
dst.delete();

最後に

関数名前にXmeansと付けていますが、最初にも言った通りXmeansは無関係なので、使うときは関数名を変えてください。