超像素经典算法SLIC的代码的深度优化和分析。
超像素经典算法SLIC的代码的深度优化和分析。
现在这个社会发展的太快，到处都充斥着各种各样的资源，各种开源的平台，如github,codeproject,pudn等等，加上一些大型的官方的开源软件，基本上能找到各个类型的代码。很多初创业的老板可能都曾经说过基本上我的程序员不需要自己写算法，但是他们要学会搜索，强有力的搜索能力基本能解决可能会遇到的一切问题，比如前一段时间流行的prisma滤镜，现在你在github上能找到一大堆了。
可我想表达的并不是这个，上述这个情况确实是正确的。但是，我相信这个世界上90%的开源代码是垃圾代码，还有9%的代码是理想代码，能够有1%的是优质代码就已经是相当不错了。基本上我们拿到的网络中的某个参考代码都还是要经过自己的细心改良和优化才能真正的应用于项目中，而这部分能力并不是每个人都有。
超像素经典的算法SLIC就属于上述1%的一员，他有论文的介绍原理性的东西，有数学公式的推导，有和其他算法的比较数据，更重要的是他还有和论文完全对应的参考代码，而且有C++、matlab以及GPU版本，可以说是非常有良心的一篇论文。虽然是优质代码，但是当你真正的去研究他的代码时，你就会发现离实际的应用距离还有很远的路要走：可怕的内存占用，大量的浮点计算还是很客观的时间开销。你在网络上搜索，包括github上，可以找到一些相关的用SLIC做图像分割的代码，而百度上所搜SLIC也能找到一些博客园或者CSDN的博客的介绍，不过大部分都停留在对源代码的解释上。
我这里不对论文的思想做详细的解释，相关的论文和源代码可以参考：http://ivrl.epfl.ch/supplementary_material/RK_SLICSuperpixels/index.html
为了描述方便，这里贴下原文对算法流程的描述：
我们来一条一条的优化和整理。
首先，论文提到图像信息是取自于CIELAB空间而不是RGB空间，在论文中有相关的描述如下：CIELAB color space is widely considered as perceptually uniform for small color space。我们实际编程测试时也发现在相同的参数下，很明显CIELAB空间的分割效果明显要比RGB空间的规则很多（如下两图所示）。这个可以从LAB空间的球形模型（左图）得到阐释吧。但是也不是说RGB空间的结果就完全不行，我们如果适当的调大M参数的值，也能获得不错的结果。唯一不同的时，使用LAB空间在时间和空间上会稍微多一点。
接着，我们来看看作者的代码，首先是RGB2XYZ、RGB2LAB、DoRGBtoLABConversion三个函数，很明显，这是RGB空间转换为LAB空间的算法，XYZ空间只是作为两者转换的中间件存在的，作者代码里采用了double类型数据来保存转换的结果，这是严格意义上的转换，作为论文的算法验证来说，是无可厚非的，但是拿到工程中使用可以说是一场灾难。8倍于原始图像数据的额外的内存占用在很多场景下就已经无法使用了，还有可怕的浮点计算的硬件屏障（不要以为大家都是在PC上使用图像^^），我们必须想办法降低这个内存占用和复杂度。
在我的博文 颜色空间系列2: RGB和CIELAB颜色空间的转换及优化算法 中，提出了一种快速算法，可以无任何浮点计算快速的将RGB转换到和原图占用内存一样大小的内存空间中，而后续的编码也证明这种转换的精度损失对于结果的影响是完全在可以接受的范围内的。
LAB空间结果
RGB空间结果
下一步就是初始化聚类中心，对应的函数在GetKValues_LABXYZ或者GetLABXYSeeds_ForGivenStepSize中，按照用户提供的超像素的大小或者超像素的个数来均匀的分布取样的XY位置和LAB的值，作者用vector来保存这些值，在源代码的其他很多地方，作者也使用了vector，就我看来，如果不是迫不得已，我基本不使用这个东西。特别作为本例，并没有使用到vector的啥高级特性，如果直接自己动态的分配内存，在很多方面都有着特别的灵活性，而且封装的东西效率很多情况下并不如意。
同样的道理，从速度方面考虑，再一次，我们也不使用源代码中的double类型来保存中心点的坐标，而是使用整形，是否可行，一切皆以最后的结果说话，实际就证明此法可行。在GetKValues_LABXYZ最后，作者为了避免初始中心点落入图像的边缘或者噪音点，进行了一次重定位操作，即取初始中心点周边3*3领域梯度值最小的那个像素为新的中心点，原文的话为：The centers are moved to seed locations corresponding to the lowest gradient position in a 3 × 3 neighborhood. This is done to avoid centering a superpixel on an edge, and to reduce the chance of seeding a superpixel with a noisy pixel 。坦白的说，个人这个其实没有啥必要，因为一般情况下图像的边缘宽度也非一个像素宽，而噪音一般也非一个孤立点，但是无论如何，做一下也好。
于是就涉及到了图像梯度值的计算了，其实说白了也就是图像的边缘算法，对应于源代码的DetectLabEdges函数，注意这里是在LAB空间的边缘检测了。哇，又是一堆的浮点计算。实际上，计算这个边缘信息只用L通道的信息就完全足够了，L就是图像的明度吗，也就表面了一个像素的强度信息。鉴于前面我们使用的整形的数据，这里我给出一个简单的单通道的边缘计算方式，还支持In-Place操作，以及对图像四周的一圈像素也是进行了计算的。
// 基于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);
}至于上述的最后的量化操作，也可以考虑使用abs(GradientH) + abs(GradientV)来实现。
下一步就是最核心的计算了，通过聚类的方式迭代计算新的聚类中心，具体过程见源代码中的PerformSuperpixelSLIC函数，我们这里提下优化问题。
迭代计算聚类中心是本算法的核心，而迭代的核心就是计算距离，因为SLIC的聚类的数据源是LAB和XY的综合体，而LAB和XY各自的取值范围不同，如果直接采用欧式距离计算，则会产生非常混乱的效果，因此，作者提出了如下的计算公式：
通过两个参数m和S来协调两种距离的比例分配。参数S表示XY空间的最大的可能值，这个通过用户输入的超像素大小可以自动计算出来（详见论文），而参数M为LAB空间的距离可能最大值，论文提出其可取的范围建议为[1,40]。
让我们来仔细看看上述公式，似乎运算量比较大， 有三次开平方，有除法还有多次乘法。 首先，由于只是距离比较，而不是基于距离数据的某种指标计量，因此，D'公式中的开平方计算是无需的。其实，在D’公式的开平方内部，Dc和Ds都有平方计算，因此，Dc和Ds计算式中的开平方计算是无需的。这已经大大的减少了计算量，接着，我们把D'开平方内部数据乘以m*m*s*s可以得到下式：
因此，我们只需要比较这个式子的值就是比较距离的相对大小了，前面已经说了LAB和XY我们都采用整形数据，聚类的中心也是整形，因此，如果S和M也是整形值，则上述计算的所有计算都是整形，而S值当然可以是整形，M值得范围[1,40]，当然也可以取整形了, So, 这里我们优化后就只有整形计算了，完全避免了浮点的计算。
但是要仔细的分析下上述式子的取值范围，如果超出了int32所能表达的范围，那在32位系统上就有所麻烦了，众所周知，在32位系统上int64数据类型的计算速度是不快的。当然如果是64位系统这个问题是不会存在的。我们经过上述的LAB颜色空间转换后，LAB各分量的范围都在[0,255]之间，因此,
这个值最大不会超过255 * 255 * 3，而S*S等于超像素的大小，超像素一般不会大于100*100把，再大也没啥意义了。同理
表示了超像素中心点和其他的距离，也一般不会大于100*100，而m为用户输入值，前面说了m最大取值40，综合所属上面的式子的最大值是不会超过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);几个细节描述下：
（1）上述代码有几句被我注释掉了，我的本意是一次性分配足够的内存，然后在分给其他的变量，这样有些代码写起来简洁些比如清零和释放，但是我测试时发现对我默认的那个测试图，被注释掉的代码会慢的比较明显（170ms和140ms），但是换一副图区别就不太明显了，后面我分析了下对于我那些测试图，其变量ActualSPAmount恰好等于1024，这个看似和对齐有很大的关系，不知道各位童鞋是怎么认为的。
（2）上述代码的中Distance数据我使用了unsigned类型的，这主要是为了后面能方便使用：
memset(Distance, 255, Width * Height * sizeof(int)); // 把距离设置为最大值
这一句代码，因为如果使用的是signed类型的，则只能用一个for循环给每个变量赋值为int类型的最大值了，memset肯定要比for快一些。
在迭代完成后，大部分工作就已经基本完成了，论文中提出大概10次迭代就足够了，其实我们实测，经过一次迭代就基本已经比较靠谱了，上述流程中那个计算residual error E的过程是不需要的，耗时有耗力，实际应用中，我觉得取迭代次数为4基本能平衡速度和精度了。
在最后，论文提出了一些后处理的过程，这主要是为了去除前面分割过程中的产生的一些比较小的分割块，相关的代码在EnforceLabelConnectivity中，这个算法的核心思想就是利用区域生长法找到图像中每个超像素的大小，然后把过于小的超像素合并到周边大的超像素中，这里有几个问题其实值得商榷：
（1）过小的超像素合并到周边哪个超像素中呢，论文是采用的是找到的最后的相邻的超像素，其实抑或是采用找到的第一个相信的超像素也好，都存在一个问题，这合理吗？我们看下面3个图：
局部原图
不进行EnforceLabelConnectivity前的结果
进行EnforceLabelConnectivity后
大家注意上述的第三张图，进行了EnforceLabelConnectivity后，红色圆圈内的结果，背景和那个项链已经练成了同一个区域，显然这个是不合理的，也是不正确的，对于后续的分割或者其他操作都是不利的，而未进行前，这两个个区域明显是分开的，所以这个操作到底要如何进行呢？？？？？？？？？？？？？？
（2）很多人可能都没有注意到进行了EnforceLabelConnectivity后，超像素的个数可能不会有任何减少，甚至还会出现增加的情况，这岂不逆天了。
EnforceLabelConnectivity不就是合并一些小的超像素到周边去吗，怎么可能增加呢，我为了这个问题也是困惑了很久很久，甚至有好几天走路吃饭都在想这个问题，有一次差点被车撞倒。也我曾怀疑作者的这个函数写的有问题，也用自己以前一个非常靠谱的区域生长来代替他的代码结果还是一样。后来有一天我突然恍然大悟，大家都没有错，问题是出在超像素分割个本身过程的，由于聚类过程的特性，并不能保证每一类在XY空间里都能连续的，比如上面中间那个图里，大量白色边界重叠在一起的区域里就有很多超像素在空间上已经分离了（其实这个时候就不能叫他为1个超像素了，而是2个或者多个），而在EnforceLabelConnectivity函数里，是基于区域生长的，也就是基于空间连续性做的判断，这样就出现了上述问题。
最后，这代码值得参考的还有各个标签之间分界线的绘制的技巧，对应的函数在DrawContoursAroundSegments里。
通过上述优化，加上其他的一些编程技巧（比如及时释放不在需要的内存），在整个函数的执行过程，大概需要3.5倍额外的图像内存，就可以完成整个过程，而作者提供的代码，少说也要20倍吧，呵呵，速度方面的，1000*1000的图，产生1600个超像素，迭代10次，大概需要300ms，而源代码大概需要3S, 10倍的提速，实际上，如果迭代4次，150ms就可以获得结果。
论文中提出，当m越大时，空间相似性越重要，超像素的结果就越紧凑，而m小时，则超像素则越靠近图像的边界，但是形状越不规则，这个我们也可以从下面的实验图片中看到：
局部原图
m = 5
m = 2
m = 40
论文中还提出了自动参数的实现，即m值无需用户输入，在迭代过程自动更新，我也实现了下，觉得这个意义不大，而且或拖慢计算和速度和增加内存消耗。 通常我们取m固定为32左右就可以了。
我做了UI界面，如下图，优化思路已经描述的相当清晰了，有能力者自可为之。
DEMO下载地址：http://files.cnblogs.com/files/Imageshop/Slic.rar
超像素有很多应用，其中一个方向就是抠图和图像分割，下一阶段我将抽空研究下这方面的东西。