针对SLIC源论文代码的分析与改进

发布于:2022-10-27 ⋅ 阅读:(379) ⋅ 点赞:(0)

1.

原:图像信息是取自于CIELAB空间。作者的代码使用RGB2XYZ、RGB2LAB、DoRGBtoLABConversion三个函数,这是RGB空间转换为LAB空间的算法。XYZ空间只是作为两者转换的中间件存在的。采用double类型数据保存转换结果。
缺点:8倍于原始图像数据的额外的内存——占内存

💛改进:采用一种 无任何浮点计算的 将RGB转换到和原图占用内存一样大小 的内存空间中。参考链接:颜色空间系列2: RGB和CIELAB颜色空间的转换及优化算法 - Imageshop - 博客园 (cnblogs.com)

2.

原:初始化聚类中心。对应的函数在GetKValues_LABXYZ或者GetLABXYSeeds_ForGivenStepSize中。原文按照用户提供的超像素大小/个数来均匀撒种子来获取XY位置和LAB的值,作者用vector来保存这些值。
💛改进:自己动态分配,具有灵活性。

3.

原:从速度方面考虑。原文使用double类型来保存中心点的坐标。
改进:我们用整形。

4.

原:在GetKValues_LABXYZ最后,作者为了避免初始中心点落入图像的边缘或者噪音点,进行了一次重定位操作,即取初始中心点周边3*3领域梯度值最小的那个像素为新的中心点。其中涉及了梯度值的计算,也就是图像的边缘算法,对应于源代码的DetectLabEdges函数,注意这里是在LAB空间的边缘检测了。然后又是一堆浮点计算。
💛改进:计算边缘信息只使用L通道就足够了。L表明一个像素的强度信息。

最后的量化操作,也可以考虑使用abs(GradientH) + abs(GradientV)来实现。

//    基于Sobel算子的边缘检测算法
//    支持In-Place操作,即Dest可以等于Src,当然处理后Src的数据就被改变了
//    四周边缘像素同样得到了处理
void GetEdgeGradient(unsigned char *Src, unsigned char *Dest, int Width, int Height)
{
    const int Channel = 1;
    unsigned char *RowCopy = (unsigned char *)malloc((Width + 2) * 3 * Channel);
    unsigned char *SqrValue = (unsigned char *)malloc(65026 * sizeof(unsigned char));
    unsigned char *First = RowCopy;
    unsigned char *Second = RowCopy + (Width + 2) * Channel;
    unsigned char *Third = RowCopy + (Width + 2) * 2 * Channel;
    
    for (int Y = 0; Y < 65026; Y++)    SqrValue[Y] = (int)(sqrtf(Y * 1.0) + 0.49999f);

    memcpy(Second, Src, Channel);
    memcpy(Second + Channel, Src, Width * Channel);                                                    //    拷贝数据到中间位置
    memcpy(Second + (Width + 1) * Channel, Src + (Width - 1) * Channel, Channel);

    memcpy(First, Second, (Width + 2) * Channel);                                                    //    第一行和第二行一样

    memcpy(Third, Src + Width * Channel, Channel);                                                    //    拷贝第二行数据
    memcpy(Third + Channel, Src + Width * Channel, Width * Channel);
    memcpy(Third + (Width + 1) * Channel, Src + Width * Channel + (Width - 1) * Channel, Channel);

    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char *LinePD = Dest + Y * Width;
        if (Y != 0)
        {
            unsigned char *Temp = First; First = Second; Second = Third; Third = Temp;
        }
        if (Y == Height - 1)
        {
            memcpy(Third, Second, (Width + 2) * Channel);
        }
        else
        {
            memcpy(Third, Src + (Y + 1) * Width * Channel, Channel);
            memcpy(Third + Channel, Src + (Y + 1) * Width * Channel, Width * Channel);                            //    由于备份了前面一行的数据,这里即使Src和Dest相同也是没有问题的
            memcpy(Third + (Width + 1) * Channel, Src + (Y + 1) * Width * Channel + (Width - 1) * Channel, Channel);
        }
        for (int X = 0; X < Width; X++)
        {
            int GradientH, GradientV;
            if (X == 0)
            {
                GradientH = First[X + 0] + First[X + 1] + First[X + 2] - (Third[X + 0] + Third[X + 1] + Third[X + 2]);
            }
            else
            {
                GradientH = GradientH - First[X - 1] + First[X + 2] + Third[X - 1] - Third[X + 2];
            }
            GradientV = First[X + 0] + Second[X + 0] * 2 + Third[X + 0] - (First[X + 2] + Second[X + 2] * 2 + Third[X + 2]);
            int Value = (GradientH * GradientH + GradientV * GradientV) >> 1;
            if (Value > 65025)
                LinePD[X] = 255;
            else
                LinePD[X] = SqrValue[Value];
        }
    }
    free(RowCopy);
    free(SqrValue);
}

5.通过聚类的方式迭代计算新的聚类中心。原:PerformSuperpixelSLIC函数。作者提出的公式:

 

 LAB颜色距离(dc)和XY空间距离(ds)两种距离是使用参数m(随图片、聚类的不同而不同,通常取1~40之间的数字)和s(超像素之间的间距)来协调分配的。计算量较大。
💛改进:只是距离比较而不是基于距离数据的某种指标计量。因此D’中的开平方无需计算。在D’公式的开平方内部,Dc和Ds都有平方计算,因此,Dc和Ds计算式中的开平方计算是无需的。——》减少了计算量——》接着可以把D’开平方内部的数据乘以m*m*s*s,可以得到:

 因此,这个式子的值就是比较距离的相对大小了。前面说了LAB和XY都采用整形数据,故聚类中心也是整形,如果S和M的值都是整形,则所有的计算就都是整形了。——》优化后所有的计算都是整形,避免了浮点型数据的计算。

注意:不能超出int32所表达的范围。在上述LAB颜色空间转换后,LAB各分量范围都在[0,255]之间,故dc的平方这个值不会超过255 * 255 * 3。S*S等于超像素的大小,一般超像素大小不会超过100*100。 同理ds的平方表示超像素中心和其他点的距离,一般不会超过100*100。综上所述,上面的式子的最大值不会超过int32所表达的最大范围。具体代码如下:

 *******  迭代更新中心点 ******** //
    
    int PowerStep = Step * Step;
    int PowerM = M * M;
    
    int *SumL = (int *)_aligned_malloc(ActualSPAmount * sizeof(int), 16);
    int *SumA = (int *)_aligned_malloc(ActualSPAmount * sizeof(int), 16);
    int *SumB = (int *)_aligned_malloc(ActualSPAmount * sizeof(int), 16);
    int *SumX = (int *)_aligned_malloc(ActualSPAmount * sizeof(int), 16);
    int *SumY = (int *)_aligned_malloc(ActualSPAmount * sizeof(int), 16);
    int *SumC = (int *)_aligned_malloc(ActualSPAmount * sizeof(int), 16);

    //int *SumMem = (int *)_aligned_malloc(ActualSPAmount * 6 * sizeof(int), 16);
    //int *SumL = SumMem + ActualSPAmount * 0, *SumA = SumMem + ActualSPAmount * 1, *SumB = SumMem + ActualSPAmount * 2;
    //int *SumX = SumMem + ActualSPAmount * 3, *SumY = SumMem + ActualSPAmount * 4, *SumC = SumMem + ActualSPAmount * 5;

    unsigned short *TempLabel = (unsigned short *)_aligned_malloc(Width * Height * sizeof(unsigned short), 16);
    unsigned int *Distance = (unsigned int *)_aligned_malloc(Width * Height * sizeof(unsigned int), 16);

    for (int Iteration = 0; Iteration < Iter; Iteration++)
    {
        memset(Distance, 255, Width * Height * sizeof(int));            //    把距离设置为最大值
        for (int Index = 0; Index < ActualSPAmount; Index++)            //    遍历所有的超像素
        {
            int CenterL = SeedL[Index];
            int CenterA = SeedA[Index];
            int CenterB = SeedB[Index];
            int CenterX = SeedX[Index];
            int CenterY = SeedY[Index];
            int X1 = CenterX - Step;                            //    以当前超像素的中心点位中心,向四周扩散Step距离
            int Y1 = CenterY - Step;
            int X2 = CenterX + Step;
            int Y2 = CenterY + Step;
            if (X1 < 0)            X1 = 0;                            //    注意不要越界
            if (X2 > Width)        X2 = Width;
            if (Y1 < 0)            Y1 = 0;

            if (Y2 > Height)    Y2 = Height;

            for (int Y = Y1; Y < Y2; Y++)
            {
                int Pos = Y * Width + X1;
                int PowerY = (Y - CenterY) * (Y - CenterY);
                for (int X = X1; X < X2; X++)
                {
                    int CurrentL = L[Pos];
                    int CurrentA = A[Pos];
                    int CurrentB = B[Pos];
                    int DistLAB = (CurrentL - CenterL) * (CurrentL - CenterL) + (CurrentA - CenterA) * (CurrentA - CenterA) + (CurrentB - CenterB) * (CurrentB - CenterB);
                    int DistXY = (X - CenterX) * (X - CenterX) + PowerY;
                    int Dist = PowerStep * DistLAB + DistXY * PowerM;        //    整形化
                    if (Dist < Distance[Pos])
                    {
                        Distance[Pos] = Dist;        //    记录最小的距离
                        TempLabel[Pos] = Index;
                    }
                    Pos++;
                }
            }
        }
        memset(SumL, 0, ActualSPAmount * sizeof(int));            //    把所有的累加数据清零
        memset(SumA, 0, ActualSPAmount * sizeof(int));
        memset(SumB, 0, ActualSPAmount * sizeof(int));
        memset(SumX, 0, ActualSPAmount * sizeof(int));
        memset(SumY, 0, ActualSPAmount * sizeof(int));
        memset(SumC, 0, ActualSPAmount * sizeof(int));

        //memset(SumMem, 0, ActualSPAmount * 6 * sizeof(int));            //    把所有的累加数据清零
        for (int Y = 0; Y < Height; Y++)
        {
            int Pos = Y * Width;
            for (int X = 0; X < Width; X++)
            {
                int Index = TempLabel[Pos];
                SumL[Index] += L[Pos];
                SumA[Index] += A[Pos];
                SumB[Index] += B[Pos];                //    累加
                SumX[Index] += X;
                SumY[Index] += Y;
                SumC[Index]++;
                Pos++;
            }
        }
        for (int Index = 0; Index < ActualSPAmount; Index++)
        {
            int Amount = SumC[Index];
            if (Amount == 0)    Amount = 1;
            SeedL[Index] = SumL[Index] / Amount;
            SeedA[Index] = SumA[Index] / Amount;    //    重新计算新的中心点
            SeedB[Index] = SumB[Index] / Amount;
            SeedX[Index] = SumX[Index] / Amount;
            SeedY[Index] = SumY[Index] / Amount;
        }
    }

    _aligned_free(L);                //    减少综合内存占用量
    _aligned_free(A);
    _aligned_free(B);
    _aligned_free(SumL);
    _aligned_free(SumA);
    _aligned_free(SumB);
    _aligned_free(SumX);
    _aligned_free(SumY);
    _aligned_free(SumC);
    //_aligned_free(TempLabel);
    _aligned_free(Distance);

    //_aligned_free(SumMem);

注意:①上面代码有部分被注释掉了,本意是一次性分配足够的内存,然后再分给其他变量。
②上述代码的中Distance数据使用了unsigned类型的,这主要是为了后面能方便使用:
 memset(Distance, 255, Width * Height * sizeof(int)); // 把距离设置为最大值

这一句代码,因为如果使用的是signed类型的,则只能用一个for循环给每个变量赋值为int类型的最大值了,memset肯定要比for快一些。

6.

论文的后处理过程——去除前面分割中产生的一些比较小的分割块,相关代码在EnforceLabelConnectivity中。这个算法的核心思想是利用区域生长法找到图像中每个超像素的大小,然后把过于小的超像素合并到周边大的超像素中。出现的问题:过小像素合并到周边哪个?论文采用的是相邻的超像素。——》会出现本是两个类别的东西分到同一类里。——》进行了EnforceLabelConnectivity后,超像素个数不会减少甚至是增加了。——》原因:问题是出在超像素分割个本身过程的,由于聚类过程的特性,并不能保证每一类在XY空间里都能连续的。比如图里,大量白色边界重叠在一起的区域里就有很多超像素在空间上已经分离了(其实这个时候就不能叫他为1个超像素了,而是2个或者多个),而在EnforceLabelConnectivity函数里,是基于区域生长的,也就是基于空间连续性做的判断,这样就出现了上述问题。

 7.DrawContoursAroundSegments里有各个标签之间分界线的绘制的技巧。
改进:及时释放不需要的内存。论文中提出,当m越大时,空间相似性越重要,超像素的结果就越紧凑,而m小时,则超像素则越靠近图像的边界,但是形状越不规则。

 8. 论文中还提出了自动参数的实现,即m值无需用户输入,在迭代过程自动更新,我也实现了下,觉得这个意义不大,而且或拖慢计算和速度和增加内存消耗。 通常我们取m固定为32左右就可以了。


SLIC源代码过程及函数

1.设定期望分割的超像素数目,打开图片。将彩色RGB图片转换为LAB空间及x、y像素坐标共5维空间。

2.DetectLabEdges。求图片中所有点的梯度=dx+dy.其中

dx=(l(x-1)-l(x+1))*(l(x-1)-l(x+1))+(a(x-1)-a(x+1))*(a(x-1)-a(x+1))+(b(x-1)-b(x+1))*(b(x-1)-b(x+1));

dy=(l(y-1)-l(y+1))*(l(y-1)-l(y+1))+(a(y-1)-a(y+1))*(a(y-1)-a(y+1))+(b(y-1)-b(y+1))*(b(y-1)-b(y+1));

3.GetLABXYSeeds_ForGivenK。给定了要分割的超像素总数K,根据LABXY信息获得种子点。

①超像素的种子点间步长Step=sqrt(N/K)。初始化种子点。按照步长均匀播撒种子点,初始化后种子点是均匀分布的。——》红色为种子点

② PerturbSeeds。扰乱种子点。在每个种子点的3*3邻域内,计算该种子点的8个邻域内像素点的Lab颜色梯度(同上述步骤2),分别与初始种子点梯度进行比较,取梯度值最小(最“平坦”)的点,并记录其LABXY信息作为新的种子点。——》绿色为更新过梯度最小值的种子点

 4.超像素的步长Step=sqrt(N/K)+2。加了一个小偏置2是为了避免Step太小,造成超像素太密集的情况。

5.PerformSuperpixelSegmentation_VariableSandM。对于每个超像素,最大的颜色距离M取值范围[1,40],一般取10。最大空间距离取步长为Step。

① 搜索范围2step* 2step(2s*2s),即设置offset=step。 在步长较短时(step<10)可以扩展offset=step*1.5作为搜索范围。

② 初始化distlab、distxy、distvec为无穷大。

maxlab初始化为10*10,maxxy初始化为step*step(s*s)

distlab代表某点与种子点的lab颜色空间距离。计算如下:

distlab(i)=(l-kseedsl(n))*(l-kseedsl(n))+(a-kseedsa(n))*(a-kseedsa(n))+(b-kseedsb(n))*(b-kseedsb(n));

distxy代表某点与种子点的空间坐标距离,计算如下:distxy(i)=(x-kseedsx(n))*(x-kseedsx(n))+(y-kseedsy(n))*(y-kseedsy(n));

dist代表某点与种子点的综合距离(归一化的颜色距离+空间距离),计算如下:dist=distlab/( maxlab)+ distxy/(maxxy);

在此提醒一下:如果将C++程序转为matlab代码时特别要注意数据类型。uint16类型变量减去double类型变量的结果是uint16类型,所以如果后者值大于前者,结果就为0。此处容易出错,需要强制类型转换。

③ 计算搜索区域内每个点离种子点的距离dist,并将搜索区域内每个点离种子点的距离保存在distvec中。因为某点可能位于多个种子点的搜索区域,所以最后保存的是离相邻种子点最近的距离,并将该像素标号为最近种子点相同的标号。同一个超像素内所有像素的标号相同。

④  计算每个新超像素内所有像素的labxy均值和坐标重心。将坐标重心作为该超像素的新种子点位置。

⑤ 上述步骤2)到4)重复迭代10次

6. EnforceLabelConnectivity。该函数主要有几个作用:保证同一个超像素都是单连通区域;去掉尺寸过小的超像素;避免单个超像素被切割的情况。

①  先计算超像素理想面积大小:SUPSZ = sz/K = N/K;

②  有两种标号向量:上一步骤中得到的旧标号向量labels(即步骤5中得到的klabels),但其存在多连通,过小超像素等问题,需要优化。新标号向量nlabels,初始化值全为-1。

③  首先选择每个超像素的起始点(左上角的第一个点),起始点判断条件:a) 按照从左到右,从上到下的“Z”型顺序查找。b)该点在新标号向量nlabels中未被标记过(值为-1)。将其坐标保存在xvec[0],yvec[0]中。

④   记录前一个相邻超像素的标号值adjlabel。判断条件:a)在步骤3中起始点的四邻域。b)在新标号向量nlabels中被标记过(标号大于0)。记录adjlabel的目的是:如果当前超像素尺寸过小,将当前超像素标号全部用adjlabel代替。即合并到前一个相邻超像素,参考下面步骤⑥。

⑤   扩展当前超像素。首先选择起始点作为当前操作的中心点,然后对其四邻域进行判断是否属于该超像素成员。判断条件:a) 该点在新标号向量nlabels中未被标记过(值为-1);b)该点n和当前操作中心点c在旧标号向量中标号一致,即labels(n)= labels(c),可以理解为原先就是属于同一个超像素的成员。如果判断是超像素的新成员,那么把该新成员作为新的操作中心点,循环直到找不到新成员为止。

⑥  如果新超像素大小小于理想超像素大小的一半(可以根据需要自己定义),将该超像素标号用前一个相邻超像素的标号值adjlabel代替,并且不递增标号值。

⑦  迭代上述步骤3)到6)直到整张图片遍历结束。

7. 绘制分割结果,退出窗口。

本文含有隐藏内容,请 开通VIP 后查看

网站公告

今日签到

点亮在社区的每一天
去签到