Edge Detection and Image SegmentatiON (EDISON) System

Edge Detection and Image SegmentatiON (EDISON) System

转载地址：http://wenku.baidu.com/view/699e0f727fd5360cba1adb46.html

Edge Detection andImage SegmentatiON (EDISON) System

 

一、概述

 

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(inti = 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维桶缓存(3dbuckets[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. 区域邻接矩阵(RegionAdjacency 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 = newRAList [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 = newint [totalBoundaryCount]; 

int*boundaryIndex = new int [regionCount];

intcounter = 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中第一个点地址

intnumRegions = 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];

}