65.9K
CodeProject 正在变化。 阅读更多。
Home

实现主成分分析图像分割

starIconstarIconstarIconstarIconstarIcon

5.00/5 (5投票s)

2024年6月21日

CPOL

3分钟阅读

viewsIcon

4709

downloadIcon

120

如何实现主成分分析聚类算法来执行图像分割。

引言

通常需要将图像分割成颜色相似的区域。这可以用作 AI 和图像识别中的一个低级步骤,但显而易见的直接应用是选择一个颜色调色板,将 24 位 rgb 真彩色图像转换为 256 色的调色板 .gif。

背景

标准算法是 k 均值聚类,这是一种用于查找给定任意距离函数的聚类数据的通用方法。主成分分析是一种新方法,它使用图像只有三个变化分量(红、绿、蓝,或在其他空间中的等效分量)的信息。方差的主成分是在 RGB(为了便于讨论)颜色空间中的轴,如果我们用第一个和最后一个数据点来切割它,它将是最长的。通过获取该分量,我们通过找到特征向量来实现,有效地拥有了一个颜色空间中立的系统。为了找到两个聚类,我们沿着垂直于该主成分的平面进行切割。找到切割位置的方法不是取平均值,而是使用 Otsu 阈值来找到最自然地划分为两种类型的位置。在实现了二分之后,执行进一步的二分就非常简单了。

使用代码

特征向量是矩阵的一个属性,例如当特征向量乘以矩阵时,特征向量保持不变,除了缩放。通常,特征向量被归一化,而缩放称为特征值。它们具有许多数学性质。但我们感兴趣的是,当矩阵是等长的一系列列表的协方差矩阵时,特征向量表示方差的主要轴,如特征值所示。因此,RGB 空间中最大的特征向量是沿其切割的向量。

由于像素只有三个分量,我们的矩阵只有 3x3,这意味着我们可以使用特殊且非常高效的代码来找到特征向量。

 

/* Eigen-decomposition for symmetric 3x3 real matrices.
   Public domain, copied from the public domain Java library JAMA. */

#ifndef _eig_h

/* Symmetric matrix A => eigenvectors in columns of V, corresponding
   eigenvalues in d. */
void eigen_decomposition(double A[3][3], double V[3][3], double d[3]);

#endif

现在我们有了一种获取特征向量的方法,我们可以对 rgb 空间中的像素进行主成分分析。

/*
   Get principal components of variance
   Params: ret - return for components of major axis of variance
           pixels - the pixels
           N - count of pixels
 */
static int pca(double *ret, unsigned char *pixels, int N)
{
    double cov[3][3];
    double mu[3];
    int i, j, k;
    double var;
    double d[3];
    double v[3][3];
    
    for(i=0;i<3;i++)
        mu[i] = meancolour(pixels, N, i);
    
    /* calculate 3x3 channel covariance matrix */
    for(i=0;i<3;i++)
        for(j=0;j<=i;j++)
        {
            var  =0;
            for(k=0;k<N;k++)
            {
                var += (pixels[k*3+i] - mu[i])*(pixels[k*3+j] - mu[j]);
            }
            cov[i][j] = var / N;
            cov[j][i] = var / N;
        }
    eigen_decomposition(cov, v, d);
    /* main component in col 3 of eigenvector matrix */
    ret[0] = v[0][2];
    ret[1] = v[1][2];
    ret[2] = v[2][2];
    
    return 0;
    
}

这应该非常简单,只需正确设置对 eigen_decomposition 的调用即可。

 

然后是核心部分。首先是 PALENTRY 结构。

    typedef struct
{
    int start;
    int N;
    double error;
    unsigned char red;
    unsigned char green;
    unsigned char blue;
} PALENTRY;

start 是像素数组的索引,N 是计数。除此之外,这里没有什么令人惊讶的。我们可以使用任何颜色距离函数来计算误差。

现在我们实际分割条目。

    /*
   split an entry using pca and Otsu thresholding.
   We find the pricipal component of variance in rgb space.
   Then we apply Itsu thresholding along that axis, and cut.
   We partition using one pass of quick sort.
 */
static void splitpca(PALENTRY *entry, PALENTRY *split, unsigned char *buff)
{
    int low, high;
    int cut;
    double comp[3];
    unsigned char temp;
    int i;
    
    pca(comp, buff + entry->start*3, entry->N);
    cut = getOtsuthreshold2(buff+entry->start*3, entry->N, comp);
    low = 0;
    high = entry->N -1;
    while(low < high)
    {
        while(low < high && project(&buff[((entry->start+low)*3)], comp) < cut)
            low++;
        while(low < high && project(&buff[((entry->start+high)*3)], comp) >= cut)
            high--;
        if(low <= high)
        {
            for(i=0;i < 3;i++)
            {
              temp = buff[(entry->start+low)*3+i];
              buff[(entry->start+low)*3+i] = buff[(entry->start+high)*3+i];
              buff[(entry->start+high)*3+i] = temp;
            }
        }
        low++;
        high--;
    }
    split->start = entry->start + low;
    split->N = entry->N -low;
    entry->N = low;
    calcerror(entry, buff);
    calcerror(split, buff);
    
}

最后是 Otsu 阈值处理。

    /*
 get the Otusu threshold for image segmentation
 Params: gray - the grayscale image
 width - image width
 height - uimage height
 Returns: threshold at which to split pixels into foreground and
 background.
 */
static int getOtsuthreshold2(unsigned char *rgb, int N, double *remap)
{
    int hist[1024] = {0};
    int wB = 0;
    int wF;
    float mB, mF;
    float sum = 0;
    float sumB = 0;
    float varBetween;
    float varMax = 0.0f;
    int answer = 0;
    int i;
    int k;
    
    for(i=0;i<N;i++)
    {
        int nc = rgb[i*3] * remap[0] + rgb[i*3+1] * remap[1] + rgb[i*3+2] * remap[2];
        hist[512+nc]++;
    }
    
    /* sum of all (for means) */
    for (k=0 ; k<1024 ; k++)
        sum += k * hist[k];
    
    for(k=0;k<1024;k++)
    {
        wB += hist[k];
        if (wB == 0)
            continue;
        
        wF = N - wB;
        if (wF == 0)
            break;
        
        sumB += (float) (k * hist[k]);
        
        mB = sumB / wB;            /* Mean Background */
        mF = (sum - sumB) / wF;    /* Mean Foreground */
        
        /* Calculate Between Class Variance */
        varBetween = (float)wB * (float)wF * (mB - mF) * (mB - mF);
        
        /* Check if new maximum found */
        if (varBetween > varMax)
        {
            varMax = varBetween;
            answer = k;
        }
        
    }
    return answer - 512;
}

我们必须稍微调整标准的 Otsu 函数,因为我们希望在特征空间中进行切割,而不是像我们通常做的那样沿着亮度进行切割。所以我们传递一个 remap 参数。 

结果

标准测试图像可从南加州大学获得。

https://sipi.usc.edu/database/database.php?volume=misc&image=11#top

 

Original 256 色
mandrill 24 bit madrill 256 color
airplane 24 bit airplane 256 colour

将结果与其他算法进行比较有点不公平。一个成熟的算法将被调整和改进许多次,并且可能产生更好的结果,尽管其核心的底层方法较弱。但是这些结果非常好,并且绝对可用。 

关注点

当显示器变为 24 位时,调色板选择似乎已经消亡。但现在它随着动画 gif 的回归以及对逼真电影的需求而复兴。

历史

 

© . All rights reserved.