Mean Shift Code for the Edge Detection and Image SegmentatiON system

【转】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];

}