【转】Mean Shift Code for the Edge Detection and Image SegmentatiON system
Edge Detection and Image SegmentatiON (EDISON) System
http://coewww.rutgers.edu/riul/research/code/EDISON/index.html
一、概述
MeanShift并不算一种很新的特征空间分析算法，但是它原理简单，计算速度较快，通常能在一次分割后形成大量小的模态区域。这样便直接将问题分析层次从像素域提升到特征域，对后续处理有很大的好处。CVPR07不少新颖的分析算法(比如多目标分割)都是以mean shift为基础的。因此，它仍然有很大的研究价值。
Rutgers的RIUL实验室将mean shift和synergistic分割算法以C++实现，并将派生的边缘检测方法集成到EDISON分析平台中，以自由软件的形式发放。本日志不讨论meanshift原理和性能，而是分析EDISON控制台程序中mean shift分割算法的实现过程和技巧。
EDISON控制台程序模块：
1. 脚本解释器(parser.h/parser.cpp/globalfnc.cpp)
由于程序参数是以脚本文件提供的，所以需要进行词法、语法分析。这不是算法的重点，这里不讨论其实现方法。调用函数为
CheckSyntax() 脚本文件语法分析，查找是否有错误语法
Run() 脚本执行
2. 算法控制平台(edison.h/edison.cpp)
控制输入输出、所有参数设置及算法执行，一般由globalfnc.cpp中EXECUTE()函数调用
3. mean shift算法(ms.h/ms.cpp/msImageProcessor.cpp)
算法核心，ms.h/ms.cpp定义了MeanShift基类，使用lattice迭代计算实现。msImageProcessor派生至MeanShift，实现了区域合并、剔除、边界查找等应用。
分割过程：
1.LoadImage 获取height, width, 数据指针pImg, 数据通道数（彩色为3，灰度为1）。
EDISON原系统仅支持PPM，PGM，PBM三种图像格式，需注意，edison不支持photoshop输出的PPM图像(ps将height width作为两行参数写入文件头；而edison默认为一行，并以空格隔开，所以需要略为修改)。我们可以很容易添加对DIB和JPG等格式的支持。
2.指定meanshift参数：
(1)spatial Bandwidth (float) 空间窗
(2)range Bandwidth (float) 特征空间窗
(3)min Region Area (int) 允许的最小区域
关于空间窗和特征空间窗的物理含义请查看Comaniciu的《Mean Shift:A Rubost Approach Toward Feature Space Analysis》。允许的最小区域对分割算法本身是无意义的，主要用于后续区域合并。
3.流程：
(1) 实例化类msImageProcessor
(2) msIP::DefineImage(pImage, channel, height, width) 定义图像数据
内部流程为
a. MS::DefineLInput() 设置lattice格形数据
b. 若用户未定义核函数，定义核函数 MS::DefineKernel()。算法默认使用单位均匀核函数(flat kernel function)
(3) msIP::Filter(spatialBW, rangeBW, speedUpLevel) 处理
其中speedUpLevel为速度等级，值为 NONE, MEDIUM, HIGH。
内部流程为：
a.根据speedUpLevel分别调用NewNonOptimizedFilter(), NewOptimizedFilter1(), NewOptimizedFilter2()进行meanshift迭代计算,NONE速度最慢，但精度最高；HIGH反之。具体区别将在以后说明。
b.Connect() 对每个像素进行模式标注
(4) msIP::GetResults(filtImage) 获取粗分割后的图像, filtImage必须预先分配内存
(5) msIP::FuseRegions(spatialBW, miniRegion) 根据设定的分割最小区域合并
它包括区域邻接矩阵RAM建立、传递闭包搜索和小区域剔除。此过程较为复杂，且不属于mean shift本身，将在后续作详细分析。
(6) msIP::GetResults(segmImage) 获取合并区域后的最终结果, segmImage必须预先分配内存
(7) RegionList *regionList = msIP::GetBoundaries() //获取以区域边界点为存储对象的区域数组
int *regionIndeces = regionList->GetRegionIndeces(0);
int numRegions = regionList->GetNumRegions(); //获取区域总数
numBoundaries_ = 0;
for(int i = 0; i < numRegions; i++)
numBoundaries_ += regionList->GetRegionCount(i); //获取边界点总数
boundaries_ = new int [numBoundaries_];
for(i = 0; i < numBoundaries_; i++)
boundaries_[i] = regionIndeces[i]; //赋予边界点在原始图像数据的索引值
(8) 存储分割后数据,数据指针为segmImage
(9) 存储粗分割数据,数据指针为filtImage
(10) 存储带边界的分割数据，须在segmImage上将所有boundaries_[i]设置为边界颜色。
二、迭代核
设图像数据长度为L，通道数为3(LUV格式存储)，数据指针为data，空间窗为sigmaS，特征空间窗为sigmaR。以无加速(NO_SPEEDUP)设置下的NewNonOptimizedFilter()说明edsion迭代原理。
非参数核方法的关键为窗的选取以及窗内点选取的代码实现。如果采用常规方法，那么需要首先提取当前点迭代P所在窗以及其邻接四个窗内的所有点，然后比较每个点和P的特征分量距离是否在特征空间窗内。这个过程中的比较次数为O(L*I*sigmaS*sigmaS)，I为每个数据点的平均迭代次数，和图像特征有关。同时，迭代过程需要频繁进行数据指针移动和边界检查，很容易造成计算错误。
edison采用了一种3维桶缓存(3d buckets[x,y,L])来替换搜索，其过程如下：
(1)分别对每个数据点在空间域(x,y)和特征空间域(L,u,v)加窗，装入数据缓存sdata
for(i=0; i<L; i++)
{
sdata[idxs++] = (i%width)/sigmaS;
sdata[idxs++] = (i/width)/sigmaS;
sdata[idxs++] = data[idxd++]/sigmaR;
sdata[idxs++] = data[idxd++]/sigmaR;
sdata[idxs++] = data[idxd++]/sigmaR;
}
(2)计算sdata在(x,y,L)3个分量下的尺度，构造3维数组作为桶缓存，并初始化
int nBuck1, nBuck2, nBuck3;
int cBuck1, cBuck2, cBuck3, cBuck;
nBuck1 = (int) (sMaxs[0] + 3); // sMaxs[0]为x分量最大值，最小值为0
nBuck2 = (int) (sMaxs[1] + 3); // sMaxs[1]为y分量最大值，最小值为0
nBuck3 = (int) (sMaxs[2] - sMins + 3); // sMaxs[2]，sMins为L分量极值
buckets = new int[nBuck1*nBuck2*nBuck3];
for(i=0; i<(nBuck1*nBuck2*nBuck3); i++)
buckets[i] = -1;
(3)根据每个点(x,y,L)值，计算其在buckets中的对应位置，然后将点坐标装入slist中
int* slist = new int[L];
int idxs = 0;
for(i=0; i<L; i++)
{
cBuck1 = (int) sdata[idxs] + 1;
cBuck2 = (int) sdata[idxs+1] + 1;
cBuck3 = (int) (sdata[idxs+2] - sMins) + 1;
cBuck = cBuck1 + nBuck1*(cBuck2 + nBuck2*cBuck3);// 3维数组坐标
slist[i] = buckets[cBuck];
buckets[cBuck] = i;
idxs += lN;
}
可以看到，此操作完成后，buckets[x,y,L]存储了所有值为(x,y,L)的点空间坐标最大值loc，若loc_next=slist[loc]不为-1，则表示了下一个值为(x,y,L)的点。因此，通过不断访问loc_next，就能找到(x,y,L)在窗sigmaS*sigmaS*sigmaR下的所有邻接点。buckets和slist构造完成后，就可以进行迭代运算了。
buckets有点在于完全避免了O(L*I*sigmaS*sigmaS)的比较运算以及其带来的大量指针移动和边界检查运算。它的代价在于nBuck1*nBuck2*nBuck3+L的内存开销。其缺陷是由于buckets是根据各分量尺度构造的，要求sigmaS和sigmaR不变，故无法推广到自适应变带宽mean shift算法中。
(4)迭代过程按照标准mean shift算法完成，可以根据文献直接分析，这里不做阐述。值得注意的是：
double hiLTr = 80.0/sigmaR;
是高L值像素的阈值，即当值大于hiLTr时，误差倍数按乘4计算
NONE，MEDIUM，HIGH的区别
(1) MEDIUM_SPEEDUP模式
设P为当前点，每次迭代后获得meanshift向量Mh，将窗口移动到点P1=P+Mh后，并不立即进行迭代。而是计算P1[x,y]与sdata[x,y]在特征空间的差值。当差值小于阈值TC_DIST_FACTOR，按下列规则快速分配模式：
a.若sdata[x,y]已分配某一模式，则将P1标记为此模式；
b.否则，记录存在从P到sdata的移动路径至pointList，继续迭代至收敛。然后将pointList上的所有点赋予同一模式。
代码如下，modeTable为模式标志数组，0表示未分配模式，且无meanshift路径经过，1表示已分配模式，2表示有meanshift路径经过。
if (diff < TC_DIST_FACTOR)
{
if (modeTable[modeCandidate_i] == 0)
{
pointList[pointCount++] = modeCandidate_i; //加入至路径链表
modeTable[modeCandidate_i] = 2; //标记有路径经过
}
else
{
for (j = 0; j < N; j++)
yk[j+2] = msRawData[modeCandidate_i*N+j];
modeTable[i] = 1; //标记已分配模式
mvAbs = -1;
break;
}
}
MEDIUM_SPEEDUP模式速度性能与TC_DIST_FACTOR设置有关。TC_DIST_FACTOR过小，性能提升不太，过大则可能造成大量计算错误。
(2) HIGH_SPEEDUP模式
HIGH_SPEEDUP在MEDIUM_SPEEDUP模式下对运算进行了进一步精简：
在迭代过程中，若P和它邻域内某点Pn的5维差值小于阈值speedThreshold，则直接将Pn加入移动路径pointList。显然，HIGH模式能够更快速的实现收敛，但也会获得大量的错数据。
if (diff < speedThreshold)
{
if(modeTable[idxd] == 0)
{
pointList[pointCount++] = idxd;
modeTable[idxd] = 2;
}
}
迭代结束后，数据存储于LUV_data数组中。edison调用Connect()函数对各个像素进行模式标注，同时完成对各个模式的计数。Connect()将调用Fill()通过简单的非递归泛洪完成标注。
三、区域合并
区域合并包括相似的邻接区域合并和小区域剔除两个过程。它们的算法核心内容均是对区域邻接矩阵进行传递闭包迭代运算。因此，这里仅分析邻接区域合并。
1. 区域邻接矩阵(Region Adjacency Matrix, RAM)建立，函数为BuildRAM()
RAM实际上为一区域链表数组:
RAList *raList = new RAList [regionCount]; //regionCount为区域总数
它的每个元素raList[i]均为区域链表，RAList数据成员如下：
int label; //本区域标识
float edgeStrength; //边界强度，须由用户定义，一般不用
int edgePixelCount; //边界点计数
RAList *next; //指向下一个邻接区域的指针
RAM建立好后，raList[i]的所有邻接区域均按label的大小顺序链接起来。所以RAM可以被看成一个稀疏矩阵的数据结构。RAM中的所有元素均分配于raPool中：
RAList *raPool = new RAList [NODE_MULTIPLE*regionCount];
NODE_MULTIPLE=10，表示单个区域的最大邻域数。
查找规则为：
a.设当前点P区域标识为i，检查P的右边节点Pr，设其标识为j，若Pr和P的标识不同。则从raPool取出两个节点raNode1，raNode2，分别赋予标识i，j。将raNode2，raNode1分别插入raList[i]和raList[j]中。插入操作将调用RAList::Insert()，与一般的链表插入机制相同。
b.检查P的下方节点Pb，与Pr类似。
2.传递闭包迭代
由函数TransitiveClosure()完成，实际上BuildRAM()也是在TransitiveClosure内调用的。
TransitiveClosure实现过程如下：
a.检查raList[i]的每个邻近区域，若它们差异小于某阈值，则将较大的标识用较小的代替。
for(i = 0; i < regionCount; i++)
{
neighbor = raList[i].next;
while(neighbor)
{
if((InWindow(i, neighbor->label)))
{
iCanEl = i;
while(raList[iCanEl].label != iCanEl)
iCanEl = raList[iCanEl].label;
neighCanEl = neighbor->label;
while(raList[neighCanEl].label != neighCanEl)
neighCanEl = raList[neighCanEl].label;
if(iCanEl < neighCanEl)
raList[neighCanEl].label = iCanEl;
else
raList[iCanEl].label = neighCanEl;
}
neighbor = neighbor->next;
}
}
InWindow(i, neighbor->label)比较两个区域是否可以合并，阈值为0.25。
iCanEl的含义是canonical element(正则节点？)，即raList[i]的i值与其label相等的节点。BuildRAM()后，raList[i]与raList[i].label是相等的。但一些raList[i]可能会在此步骤中被赋予邻域的标识labeln(labeln一定比raList[i].label小)，以后的邻接区域标识必须要保证用最小的同类标识labeln，而不是raList[i].label来代替。
这步也就是所谓传递闭包运算的意义：如果把所有区域看作总集S，InWindow(i, neighbor->label)看成定义在S->S上的关系运算，传递闭包运算即是通过raList[iCanEl].label != iCanEl找出S中所有满足同类区域。
b.标识最小化
for(i = 0; i < regionCount; i++)
{
iCanEl = i;
while(raList[iCanEl].label != iCanEl)
iCanEl = raList[iCanEl].label;
raList[i].label = iCanEl;
}
a步骤已经完全可以保证表示最小，此步没有必要
c.区域重新计数，并重新计算模式
这步属于比较简单的计数和均值运算，不做分析。
至此便基本完成了对图像的分割运算，我们可以通过后续处理进行模式间的边界标注。
四、边界标注
GetBoundaries()函数可以获取图像边界，它将调用DefineBoundaries()定义边界。边界点将存储在boundaryMap数组中：
boundaryMap = new int [L];
若boudaryMap[i]!=-1，则i为边界，且boudaryMap[i]为i所在区域的标识。
boundaryCount = new int [regionCount];
boundaryCount[i]为第i个区域的边界点数；totalBoundaryCount为总边界点计数。
设labels为记录每个点所在区域标识的数组，则边界点按如下方式定义：
a.图像第一行和第一列所有点被标为边界，
b.若i点与它某个四邻域点区域标识不同，则记为边界,
dataPoint = i*width+j;
label = labels[dataPoint];
if( (label != labels[dataPoint-1]) || (label != labels[dataPoint+1])||
(label != labels[dataPoint-width]) || (label != labels[dataPoint+width]))
{
boundaryMap[dataPoint] = label = labels[dataPoint];
boundaryCount[label]++;
totalBoundaryCount++;
}
c.将图像最后一行和最后一列记为边界。
boundaryMap含有大量无用数据，下一步将用更紧凑的数据结构存储边界：
int *boundaryBuffer = new int [totalBoundaryCount];
int *boundaryIndex = new int [regionCount];
int counter = 0;
for(i = 0; i < regionCount; i++)
{
boundaryIndex[i] = counter;
counter += boundaryCount[i];
}
for(i = 0; i < L; i++)
{
if((label = boundaryMap[i]) >= 0)
{
boundaryBuffer[boundaryIndex[label]] = i;
boundaryIndex[label]++;
}
}
与boundaryMap相比，boundaryBuffer能够连续存储边界点，各区域边界在数组中的地址偏移量由boundaryIndex[label]指定。
DefineBoundaries()返回的区域链表regionList在最后生成，它将调用RegionList::AddRegion()函数。AddRegion的3个形参的含义分别为：
int label //区域标识
int pointCount //区域边界点总数
int *indeces //边界点在boundaryBuffer中的地址偏移量
添加方法如下：
counter = 0;
for(i = 0; i < regionCount; i++)
{
regionList->AddRegion(i, boundaryCount[i], &boundaryBuffer[counter]);
counter += boundaryCount[i];
}
AddRegion将按如下方式修改regionList：
regionList[i].label = label;
regionList[i].pointCount = pointCount;
regionList[i].region = freeBlockLoc;
for(int k = 0; k < pointCount; k++)
indexTable[freeBlockLoc+i] = indeces[i];
freeBlockLoc += pointCount;
indexTable为边界点位置索引表，freeBlockLoc指示了表中第一个未使用索引。
有了boundaryBuffer和boundaryCount之后，边界的标注就简单了：
int *regionIndeces = regionList->GetRegionIndeces(0);
//获取indexTable中第一个点地址
int numRegions = regionList->GetNumRegions(); //获取区域数
numBoundaries_ = 0;
for(int i = 0; i < numRegions; i++)
{
//获取边界点总数
numBoundaries_ += regionList->GetRegionCount(i);
}
boundaries_ = new int [numBoundaries_];
for(i = 0; i < numBoundaries_; i++)
{
//获取所有边界点索引
boundaries_[i] = regionIndeces[i];
}
