算法执行效果
相关参考资料:看着玩的。
BLOB算法简述 https://blog.youkuaiyun.com/icyrat/article/details/6594574
话说这老哥写的也太"简"了吧,完全口水话,把blob算法说的很神秘,说什么把blob算法做好了,“我们在老外面前就有的吹了”。雷刚萨斯的图像处理早在70年代就出第一版了,已经03年了,过去30年了,老外不至于一个工业应用也没有吧。
【Halcon视频教程】Blob分析的基本概念和实现流程: https://zhuanlan.zhihu.com/p/75473787 就是一个卖课的。
张冬娟, 唐万有. 基于Blob算法的烫印缺陷在线检测的研究[J]. 包装工程, 2013, 34(17):4.:
相关硬核参考资料:干货,核心
Rosenfeld, Azriel, and John L. Pfaltz. “Sequential operations in digital picture processing.” Journal of the ACM (JACM) 13.4 (1966): 471-494.
不得不说,图像处理这一块,老外起步还是早,这篇1966年的文章里面,就已经把图像处理中的blob分析和Skeleton提取讲清楚了。
Trein, J., A. Th Schwarzbacher, and B. Hoppe. “FPGA implementation of a single pass real-time blob analysis using run length encoding.” MPC-Workshop, February. 2008.
在FPGA中实现单次扫描实时Blob算法。
Zhang E H , Feng J . Run-List Based Connected Components Labeling for Blob Analysis[J]. Journal of Applied Sciences, 2008.
https://cloud.tencent.com/developer/article/1081206,本文其中一个算法就是使用这个原理实现的。
基础知识:
图像:
图像,在放大后其实是由一个个像素(下图中的格子)组成的。对于灰度图像,每个格子里面是一个数,从0到255,在内存中用一个字节描述,有些处理中会归一化到0到1的浮点数。对于彩色图像,则每个格子中有3个数,分别代表RGB通道,可以表示各种颜色。
机器视觉处理的就是这种矩阵中的数字,然后经过一系列计算,得到感兴趣的信息。
区域:
当我们想要描述图像上一个区域时,采用什么样的数据结构呢?有一种方法叫游程法,也叫游程编码(RLE) 。
描述方法为:(x1,x2,y),x采用左闭右开区间。
即:
(1,6,1),(2,5,2)
采用2个三元组即描述了如上的区域。
对应到编码实现:
struct Run
{
int nStart;
int nEnd;
int nY;
Run(int _nStart, int _nEnd, int _nY) : nStart(_nStart), nEnd(_nEnd), nY(_nY) {};
};
连通判断:4邻域和8邻域。
判断两个像素是否属于同一个区域,有两种标准,分别是4邻域和8邻域。它们的区别如下:
4邻域代表像素间距离小于等于1才认为是连通的。
8邻域代表像素间距离小于等于根号2就认为是连通的。
显然,8邻域会将上图4个像素认为是连通的,而4邻域则认为这是4个区域。
blob算法在机器视觉检测领域的地位:
机器视觉检测的一般步骤:实物-》成像-》图像-》感兴趣缺陷特征提取-》缺陷能量图-》blob算法-》缺陷blob-》深度学习或者特征的分类分级-》缺陷类型与等级-》缺陷筛选/过滤-》本地检测结果-》上传-》云端报表。
成像部分包含内容很多,涉及机械,电气,光学,自动化,仪器等专业。总的目标是将实物上感兴趣的缺陷特征呈现到图像上。
特征提取部分,主要涉及图像处理算法,传统的图像对位,阈值分割,图像比对等,还有深度学习。当前深度学习其实是降低了视觉检测行业的门槛,一个开源的模型,拿点自己的数据训练,谁都可以说我能做视觉检测,唯一的问题的算法效率太低。传统的方法,门槛也没看上去那么高,直接使用Halcon或者OpenCV等视觉库,基本的检测也能满足个七七八八,笔者在学生时代也是使用OpenCV完成了多个视觉检测与测量的项目。所以当下在工业视觉检测这一行,对算法工程师要求也不高,除非检测算法都从底层写起,但这样的公司国内一只手数的过来。
blob算法:本文的主角,其实在上面的流程图中,没有必要单独拿出来说。它的目的是从缺陷能量图上,将缺陷作为一个一个的blob提取出来。所以它的输入是一张灰度图或者二值图,输出是一堆区域(Region)。
缺陷分类分级:每个区域作为一个缺陷,具体属于什么缺陷类型,是否严重,需要使用基于特征或者深度学习的分类分级方法。
缺陷过滤与筛选,缺陷展示,缺陷上传:这部分涉及到软件。根据具体的展示需求,交互需求,通信需求等有不同的实现。
blob算法类型:
基于灰度突变,梯度大于阈值的地方作为边界,将图像分成多个blob。基于灰度突变,无法区分背景和前景,且在实际中由于阈值的选取,很难有完全闭合的边界,所以在工程上应用较少。
基于单阈值,将高于阈值的连通区域作为一个blob。本文主要介绍这种类型,有深度优先遍历和广度优先遍历两种实现方式。还会包含一些算法效率优化的技巧。
基于多阈值:局部极值的分水岭算法,使用多个阈值分割图像,将多分割图像中重心接近的blob合并。该方式比较耗时,在opencv中的实现为SimpleBlobDetector。
OpenCV的实现方式:SimpleBlobDetector
多阈值局部极值分水岭算法原理:
-
Convert the source image to binary images by applying thresholding with several thresholds from
minThreshold (inclusive) to maxThreshold (exclusive) with distance thresholdStep between
neighboring thresholds. -
Extract connected components from every binary image by findContours and calculate their
centers. -
Group centers from several binary images by their coordinates. Close centers form one group that
corresponds to one blob, which is controlled by the minDistBetweenBlobs parameter. -
From the groups, estimate final centers of blobs and their radiuses and return as locations and
sizes of keypoints.这个只是opencv中的注释,想要了解更多细节,可以参考opencv源码。
halcon的实现方法:connection
read_image (Image, 'particle')
threshold (Image, BrightPixels, 120, 255)
connection (BrightPixels, Particles)
area_center (Particles, Area, Row, Column)
通过connection函数分割区域后,通常需要后处理,比如:抑制小区域、给定方向的区域或靠近其他区域的区域。 halcon也提供了一些了区域后处理算子,比如:形态算子opening_circle 和opening_rectangle1 常用于抑制噪声,而Closing_circle 和closed_rectangle1 常用于填补空白。select_shape、select_shape_std 和 select_shape_proto 选择具有特定特征的 Blob。
可以看到,关于这个算子的原理,halcon手册上也没怎么写,我估计就是我将要介绍到的下面的两种。
深度优先的方式
算法描述
原理比较简单,直接看代码吧。
编码实现
游程定义
//游程定义,X都是左闭右开
struct Run
{
int nX1;
int nX2;
int nY;
};
//游程列表
typedef list<Run> RunList;
//每行一个游程列表
RunList* runLists = new RunList[nH];
提取游程
Mat im = imread("PDTestOut.bmp",-1);
int nW = im.cols;
int nH = im.rows;
uint8_t* pSrc = im.data;
int y, x;
for (y = 0; y < nH; y++)
{
bool bIn = false;
auto pSrc1 = pSrc + y * nW;
Run run;
for (x = 0; x < nW; x++)
{
if ((uint8_t)pSrc1[x] > 0)
{
if (!bIn) //当前不在游程内但高于阈值,记录起点
{
run.nX1 = x;
bIn = true;
}
}
else
{
if (bIn) //当前在游程但低于阈值,记录终点,push到对应行的list
{
run.nX2 = x;
run.nY = y;
runLists[y].push_back(run);
bIn = false;
}
}
}
if (bIn)
{
run.nX2 = x;
run.nY = y;
runLists[y].push_back(run);
bIn = false;
}
}
区域定义
struct Blob
{
//list<Run> runs;
Run runs[32768];
int nRunCount = 0;
int nArea;
};
typedef list<Blob> Blobs;
Blobs blobs;
游程连接
int s = 0;
bool bFinish = false;
while (!bFinish)
{
bFinish = true;
for (y = 0; y < nH; y++)
{
if (runLists[y].size() > 0)
{
bFinish = false;
//创建一个新的blob
Blob pblob;
auto& run = runLists[y].front();
pblob.runs[pblob.nRunCount++] = run;
runLists[y].erase(runLists[y].begin(), ++runLists[y].begin());
//将与该blob相连的游程全部连接
for (s = 0; s < pblob.nRunCount; s++)
{
//当前遍历到的游程
int nX1 = pblob.runs[s].nX1;
int nX2 = pblob.runs[s].nX2;
int nY = pblob.runs[s].nY;
//和上一行游程连接
if (nY > 0 && runLists[nY - 1].size() > 0)
{
for (auto it = runLists[nY - 1].begin(); it != runLists[nY - 1].end();)
{
int _nX1 = it->nX1;
int _nX2 = it->nX2;
int _nY = it->nY;
//游程都是左闭右开,8邻域
//if (_nX2 >= nX1 && _nX1 <= nX2)
//4邻域
if (_nX2 >= nX1 + 1 && _nX1 <= nX2 - 1)
{
pblob.runs[pblob.nRunCount].nX1 = _nX1;
pblob.runs[pblob.nRunCount].nX2 = _nX2;
pblob.runs[pblob.nRunCount++].nY = _nY;
runLists[nY - 1].erase(it++);
}
else
{
++it;
}
if (_nX1 >= nX2) break;
}
}
//和下一行游程连接
if (nY < nH - 1 && runLists[nY + 1].size() > 0)
{
for (auto it = runLists[nY + 1].begin(); it != runLists[nY + 1].end(); )
{
int _nX1 = it->nX1;
int _nX2 = it->nX2;
int _nY = it->nY;
//游程都是左闭右开,8邻域
//if (_nX2 >= nX1 && _nX1 <= nX2)
//4邻域
if (_nX2 >= nX1 + 1 && _nX1 <= nX2 - 1)
{
pblob.runs[pblob.nRunCount].nX1 = _nX1;
pblob.runs[pblob.nRunCount].nX2 = _nX2;
pblob.runs[pblob.nRunCount++].nY = _nY;
runLists[nY + 1].erase(it++);
}
else
{
++it;
}
if (_nX1 >= nX2) break;
}
}
}
blobs.push_back(pblob);
}
}
}
算法执行效果:
对于下图
共找到100个blob,其中第一个blob有5个游程。
给像素着色,便于可视化:
void setValue3(Mat im, Scalar scalar, Rect rc)
{
for (int y = rc.tl().y; y < rc.br().y; y++)
{
for (int x = rc.tl().x; x < rc.br().x; x++)
{
//im.at<char>(y, x) = value;
im.at<Vec3b>(y, x)[0] = scalar.val[0];
im.at<Vec3b>(y, x)[1] = scalar.val[1];
im.at<Vec3b>(y, x)[2] = scalar.val[2];
}
}
}
//给像素着色
Mat imShow(nH, nW, CV_8UC3);
imShow.setTo(0);
for (const auto& blob : blobs)
{
int nIndex = 0;
//随机生成一种颜色
int r, g, b;
do
{
r = rand() % 255;
g = rand() % 255;
b = rand() % 255;
} while (r + g + b < 200);
for (nIndex = 0; nIndex < blob.nRunCount; nIndex++)
{
int nw = blob.runs[nIndex].nX2 - blob.runs[nIndex].nX1;
Rect rc(blob.runs[nIndex].nX1, blob.runs[nIndex].nY,nw,1);
setValue3(imShow, Scalar(r, g, b), rc);
}
}
广度优先的方式
参考资料:https://cloud.tencent.com/developer/article/1081206
概述
按照处理对象的不同, 目前典型的连通性分析算法包括基于像素的方法和基于游程的方法。后者是对像素法的一种改进,它充分利用了区域各部分之间的连通关系,搜索空间得到压缩,整体效率高于前者, 因此近年来得到了更多的关注。在具体实现上,这两类方法都可采用递归法或序贯算法。递归法实现起来简单,但运行时需要消耗大量堆栈, 除了效率低,在实际应用中还容易因堆栈资源耗尽而造成算法不稳定。序贯法在扫描过程中会出现标记冲突现象,为此,常规的做法是对图像( 或子图像) 进行二次或多次扫描, 并利用冲突等价表等辅助措施来消除标记冗余 。由于等价表结构复杂,增加内存消耗; 利用等价表进行标记合并及重复扫描时间开销很大,不利于实时性要求较高的应用。
简单来说,本文采用步进式动态扫描方式,每个游程仅需扫描一次,且不必与相邻行的所有游程进行比较, 算法的搜索空间得到压缩; 游程连通性比较的分支少,简化了判断过程, 提高了操作效率; 所设计的游程及目标对象的数据结构允许由任一游程节点快速访问其所属链表的首部和尾部, 不仅为后续的数据访问提供了便利, 且提高了标记冲突时链表合并的操作速度,避免了冲突等价表的介入。实验结果表明该算法具有鲁棒、 高效的特性。
算法描述
2.1 游程及 Blob 目标对象数据结构定义
不失一般性,设分割得到的二值图像中,背景像素灰度为0,目标像素灰度为 1。一行中灰度值连续为 1 的像素构成一个游程数据单元。定义如下的游程数据结构对之进行描述:
struct RLE { short r,s,e; RLE * pn; BLOB **ppB}
其中: r 代表游程所在行号,s 为游程像素在该行的起始位置,e 为其在该行的终止位置。因每个游程数据单元必属于且仅属于某个唯一的 Blob 对象,将同属一个目标对象的所有游程数据单元组织成一个线性链表, 每个游程数据单元即为链表中的一个节点,用指针 pn 来指向链表中的下一游程节点。数据域 ppB 为间接指向其所属 Blob 对象的二级指针, 通过它可快速检索到其所属链表的头部和尾部, 用以解决标记冲突时对链表的快速合并操作, 具体过程将在后文详述。Blob 目标对象定义为:
struct BLOB { RLE * ph; RLE * pt; }
其中: ph 指向链表的头部( header) ,pt 指向链表的尾部( tail) 。可见,一个 BLOB 对象实际上描述了一个 RLE 链表, 通过它可访问同属该目标的所有 RLE 对象。算法结束后, 将动态生成一个 BLOB 链表,它描述了一幅图像中的全部目标对象。
2.2 数据准备
顺序扫描二值图像的每一行,可得到整幅图像的 RLE 表达形式。为了能够快速访问各行的游程数据, 为图像的每行维护一个一维的动态数组,数组元素类型为 RLE* ,即该行游程数据单元指针构成的索引; 若某行不存在游程数据( 即全部为背景像素) ,则数组为空。
2.3 连通性判据
相邻两行的任意两个游程连通, 当且仅当其中一个游程存在至少一个像素与另一个游程中的像素连通。游程连通性有 4 连通和 8 连通之分, 本文仅考虑 8 连通性。对于第 i - 1行和第 i 行的两个游程 RLE[ i - 1] 与 RLE[i],可用式( 1) ~( 2) 进行连通性判断若同时满足式( 1) ~ ( 2) , 则两者连通; 否则不连通[ 2] 。
RLE[ i].s≤RLE[ i - 1].e + 1 ( 1)
RLE[ i].e + 1≥RLE[ i - 1].s ( 2)
2.4 算法流程
算法约定: 设二值图像高度为 H, 记第 i 行的游程个数为size( i) 。从第 0 行开始, 按照从左到右的顺序扫描该行的游程数组。算法每次取出当前行( i) 的第 k 个及上一行( i - 1)
的第 k’个游程数据( 记为 RLE( k) , RLE( k’) , 分别称为当前游程和参考游程) 进行比较。为方便描述, 约定以符号 & 表示取对象的地址,* 表示对该地址所指向的对象进行操作。
算法完整步骤描述如下:
第 1 步 :如果 i≥H, 即图像所有行已经分析完毕, 则算法结束; 否则初始化当前游程和参考游程的索引 k、 k’←0, 转第 2 步。
第 2 步 :如果 k≥size( i) , 即当前行的游程已经分析完毕,则将行号 i 加 1,转第 1 步; 否则转第 3 步。
第 3 步 :如果 k’ ≥size( i - 1) , 说明上一行的游程已经比较完毕,则执行第 3.1 步; 否则转第 4 步。
第 3.1 步 :如果当前游程的 ppB 不为空, 说明该游程已经标记,转第 3. 2 步; 否则, 应向 BLOB 链表添加一个新的BLOB 对象及其索引 Ref←&BLOB, 并设置其 ph 和 pt 指针同时指向当前游程: BLOB.ph、 BLOB.pt←& RLE( k) ; 同时设置当前游程指向该 BLOB, 即 RLE( k) .ppB←&Ref; 继续执行第3.2 步。
第 3.2 步 :将当前游程索引值 k 增 1,转第 2 步。
第 4 步 :此时当前游程和参考游程均有效, 利用连通性判据式( 1) 和( 2) 进行两者的比较,根据比较结果按以下 3 种情况进行处理:
情况 1
如果两个游程不连通,且当前游程像素起始位置在 参 考 游 程 像 素 终 止 位 置 右 边, 即 RLE ( k ).s >RLE( k’) .e + 1,如图 1 所示,则将 k’加 1,转第 2 步。
情况 2
如果两个游程不连通, 且当前游程像素终止位置 在 参 考 游 程 像 素 起 始 位 置 左 边,即 RLE( k) .e + 1 <RLE( k’).s,如图 2 所示,则转第 3.1 步。
情况 3
如果两个游程连通,根据当前游程的 ppB 域是否为空进行以下处理,然后继续第 5 步:
ppB是否为空情况1 如果当前游程的 ppB 为空, 即其尚未标记过, 此时应直接将其挂接到上行参考游程所在链表的尾部,即:
( * RLE( k’) .ppB).pt.pn←& RLE( k)
并更新该链表的尾部,使其指向当前游程:
( * RLE( k’).ppB).pt ←& RLE( k)
同时修改当前游程所属的 ppB 指针, 使之与参考游程相同:
RLE( k).ppB ← RLE( k’).ppB
ppB是否为空情况2 如果当前游程 ppB 不为空, 即其曾经标记过, 此时应该判断当前游程所属链表是否与参考游程所属链表一致( 通过比较二者的 ppB 域是否具有相同值) 。如一致, 无需进行任何操作; 否则意味着出现标记冲突, 应合并两个链表, 为此执行下列操作:
a) 将当前游程所在的链表挂接到参考游程所在链表的尾部:
( * RLE( k’).ppB).pt.pn ←( * RLE( k).ppB).ph
b) 参考游程所在链表的尾部修改为指向当前游程所在链表的尾部:
( * RLE( k’).ppB).pt ←( * RLE( k).ppB).pt
c) 遍历 BLOB 索引数组, 将所有指向当前游程合并前所属 BLOB 的索引值修改为指向参考游程所属的 BLOB; 同时从BLOB 链表中删除当前游程在合并前所属的 BLOB 节点。
第 5 步: 更新 k 和 k’的值: 若上行参考游程像素终止位置不在 当 前 游 程 像 素 终 止 位 置 右 边, 即 RLE ( k ).e ≥RLE( k’).e,则将 k 加 1; 若参考游程像素终止位置不在当前
游程像素终止位置左边,即 RLE( k).e≤RLE( k’).e, 则将 k’加 1。返回到第 2 步继续。
算法编码实现:
数据结构定义
//数据结构定义,X都是左闭右开
struct BLOB;
struct RLE
{
int nX1 = -1;
int nX2 = -1;
int nY = -1;
RLE* pNext = nullptr; //指向区域游程链表中的下一游程节点
BLOB* ppBlob = nullptr; //指向其所属 Blob 对象的指针
};
struct BLOB
{
RLE* pHead = nullptr; //指向区域游程链表头部
RLE* pTail = nullptr; //指向区域游程链表的尾部
BLOB* pNext = nullptr;
bool isDelete = false;
};
数据准备
Mat im = imread("blob/0/100.bmp", -1);
//Mat im = imread("PDTestOut.bmp", -1);
int nW = im.cols;
int nH = im.rows;
uint8_t* pSrc = im.data;
游程提取
//提取游程
RLE** rle = new RLE*[nH];
for (int i = 0; i < nH; i++)
rle[i] = new RLE[256]; //一行最多256个游程
uint16_t* nCount = new uint16_t[nH]; //每行有多少个游程
memset(nCount, 0, nH * sizeof(uint16_t));
int y, x;
for (y = 0; y < nH; y++)
{
bool bIn = false;
auto pSrc1 = pSrc + y * nW;
for (x = 0; x < nW; x++)
{
if ((uint8_t)pSrc1[x] > 0)
{
if (!bIn) //当前不在游程内但高于阈值,记录起点
{
rle[y][nCount[y]].nX1 = x;
bIn = true;
}
}
else
{
if (bIn) //当前在游程但低于阈值,记录终点
{
rle[y][nCount[y]].nX2 = x;
rle[y][nCount[y]].nY = y;
nCount[y]++;
if (nCount[y] >= 256)
break;
bIn = false;
}
}
}
if (bIn)
{
rle[y][nCount[y]].nX2 = x;
rle[y][nCount[y]].nY = y;
nCount[y]++;
bIn = false;
}
}
游程连接
//遍历游程并连接
BLOB* pResult = nullptr; //要一直指向头部
BLOB* pBlob = nullptr; //要一直指向当前最新
for (int y = 0; y < nH; y++)
{
//Step1:初始化索引为0
int kcur = 0;
int kpre = 0;
while (1)
{
//Step2:判断当前行游程是否分析完毕
if (kcur >= nCount[y]) //当前行游程分析完毕
break;
//Step3:判断上一行游程是否分析完毕
if ((y>0 && kpre >= nCount[y - 1]) || y==0) //上一行分析完毕 或者没有上一行
{
//Step3.1:如果当前游程的 ppB 为空,应向 BLOB 链表添加一个新的BLOB 对象及其索引 Ref←&BLOB,然后转到第3.2步
if (rle[y][kcur].ppBlob == nullptr)
{
//向 BLOB 链表添加一个新的BLOB 对象及其索引 Ref←&BLOB,
if (pBlob == nullptr)
{
pResult = pBlob = new BLOB;
}
else
{
pBlob->pNext = new BLOB;
pBlob = pBlob->pNext;
}
//并设置其 ph 和 pt 指针同时指向当前游程: BLOB.ph、 BLOB.pt←& RLE( k) ;
pBlob->pHead = &rle[y][kcur];
pBlob->pTail = &rle[y][kcur];
//同时设置当前游程指向该 BLOB, 即 RLE( k) .ppB←&Ref;
rle[y][kcur].ppBlob = pBlob;
}
//Step3.2:将当前游程索引值 k 增 1,转第 2 步。
kcur++;
continue; //转第 2 步
}
//Step4:需要比较两个游程
//情况1:当前游程像素起始位置在参考游程像素终止位置右边
//RLE(k).s >RLE(k') .e + 1
//需要将上一行游程向后移动一格,继续比较
if (rle[y][kcur].nX1 > rle[y - 1][kpre].nX2)
{
kpre++;
continue; //转第 2 步
}
//情况2:当前游程像素终止位置 在参考游程像素起始位置左边
//RLE( k) .e + 1 <RLE( k').s
//需要新建一个blob
if (rle[y][kcur].nX2 < rle[y - 1][kpre].nX1)
{
//Step3.1:如果当前游程的 ppB 不为空, 说明该游程已经标记,转第 3.2 步
if (rle[y][kcur].ppBlob == nullptr)
{
//否则,应向 BLOB 链表添加一个新的BLOB 对象及其索引 Ref←&BLOB,
if (pBlob == nullptr)
{
pResult = pBlob = new BLOB;
}
else
{
pBlob->pNext = new BLOB;
pBlob = pBlob->pNext;
}
//并设置其 ph 和 pt 指针同时指向当前游程: BLOB.ph、 BLOB.pt←& RLE( k) ;
pBlob->pHead = &rle[y][kcur];
pBlob->pTail = &rle[y][kcur];
//同时设置当前游程指向该 BLOB, 即 RLE( k) .ppB←&Ref;
rle[y][kcur].ppBlob = pBlob;
}
//Step3.2:将当前游程索引值 k 增 1,转第 2 步。
kcur++;
continue;
}
//情况3:如果两个游程连通
//RLE[i].s≤RLE[i - 1].e + 1(1)
//RLE[i].e + 1≥RLE[i - 1].s(2)
//8邻域,这里我是左闭右开,所以和公式不一样
if (rle[y][kcur].nX1 <= rle[y - 1][kpre].nX2 &&
rle[y][kcur].nX2 >= rle[y - 1][kpre].nX1)
{
if (y == 185)
cout << "debug" << endl;
//根据当前游程的 ppB 域是否为空进行以下处理,然后继续第5步
if (rle[y][kcur].ppBlob == nullptr)
{
//如果当前游程的 ppB 为空, 即其尚未标记过, 此时应直接将其挂接到上行参考游程所在链表的尾部,即 :
//(*RLE(k') .ppB).pt.pn←& RLE( k)
(rle[y - 1][kpre].ppBlob)->pTail->pNext = &rle[y][kcur];
//并更新该链表的尾部,使其指向当前游程:
//(*RLE(k').ppB).pt ←& RLE( k)
(rle[y - 1][kpre].ppBlob)->pTail = &rle[y][kcur];
//同时修改当前游程所属的 ppB 指针, 使之与参考游程相同 :
//RLE(k).ppB ← RLE(k').ppB
rle[y][kcur].ppBlob = rle[y - 1][kpre].ppBlob;
}
else
{
//如果当前游程 ppB 不为空, 即其曾经标记过,
//此时应该判断当前游程所属链表是否与参考游程所属链表一致(通过比较二者的 ppB 域是否具有相同值) 。
//如一致, 无需进行任何操作; 否则意味着出现标记冲突, 应合并两个链表, 为此执行下列操作 :
if (rle[y][kcur].ppBlob != rle[y - 1][kpre].ppBlob)
{
//意味着出现标记冲突, 应合并两个链表, 为此执行下列操作 :
//a) 将当前游程所在的链表挂接到参考游程所在链表的尾部:
//(*RLE(k').ppB).pt.pn ←( * RLE( k).ppB).ph
(rle[y - 1][kpre].ppBlob)->pTail->pNext = (rle[y][kcur].ppBlob)->pHead;
//b) 参考游程所在链表的尾部修改为指向当前游程所在链表的尾部:
//(*RLE(k').ppB).pt ←( * RLE( k).ppB).pt
(rle[y - 1][kpre].ppBlob)->pTail = (rle[y][kcur].ppBlob)->pTail;
//c) 遍历 BLOB 索引数组, 将所有指向当前游程合并前所属 BLOB 的索引值修改为指向参考游程所属的 BLOB;
//同时从BLOB 链表中删除当前游程在合并前所属的 BLOB 节点。
auto pTmpBlob = pResult;
//如果是第一个,则将结果指针移动到下一个
if (pTmpBlob == rle[y][kcur].ppBlob)
{
//改变游程的指向
auto pRun = pTmpBlob->pHead;
while (pRun)
{
pRun->ppBlob = rle[y-1][kpre].ppBlob;
pRun = pRun->pNext;
}
//将pBlob首节点删除
pTmpBlob->isDelete = true;
}
else
{
//否则,需要
while (pTmpBlob)
{
if (pTmpBlob->pNext == rle[y][kcur].ppBlob)
{
//将pTmp节点删除
auto pTmp = pTmpBlob->pNext;
//改变游程的指向
auto pRun = pTmp->pHead;
while (pRun)
{
pRun->ppBlob = rle[y-1][kpre].ppBlob;
pRun = pRun->pNext;
}
pTmpBlob->pNext->isDelete = true;
break;
}
pTmpBlob = pTmpBlob->pNext;
}
}
}
}
//Step5:更新 k 和 k'的值:
//若上行参考游程像素终止位置不在 当 前 游 程 像 素 终 止 位 置 右 边, 即 RLE ( k ).e ≥RLE( k').e,
//则将 k 加 1;
int nRawkpre = kpre;
if (rle[y][kcur].nX2 >= rle[y - 1][kpre].nX2)
kpre++;
//若参考游程像素终止位置不在当前游程像素终止位置左边,即 RLE(k).e≤RLE(k').e,
//则将 k'加 1。
if (rle[y][kcur].nX2 <= rle[y - 1][nRawkpre].nX2)
kcur++;
//返回到第 2 步继续。
continue;
}
}//while 等待一行处理完成
}//行遍历
可视化展示
//可视化展示
Mat imShow(nH, nW, CV_8UC3);
imShow.setTo(0);
int nBlobCount = 0;
while (pResult)
{
if (pResult->isDelete)
{
pResult = pResult->pNext;
continue;
}
//cout << "Blob:" << nBlobCount++ << endl;
RLE* pRun = pResult->pHead;
int r, g, b;
do
{
r = rand() % 255;
g = rand() % 255;
b = rand() % 255;
} while (r + g + b < 200);
while (pRun)
{
//cout << "S:" << pRun->nX1 << ",E:" << pRun->nX2 << ",Y:" << pRun->nY << endl;
int nw = pRun->nX2 - pRun->nX1;
Rect rc(pRun->nX1, pRun->nY, nw, 1);
setValue3(imShow, Scalar(r, g, b), rc);
pRun = pRun->pNext;
}
//cout << endl;
pResult = pResult->pNext;
}
环境清理
for (int i = 0; i < nH; i++)
delete[] rle[i];
delete[] rle;
delete[] nCount;
//blob链表清理
//...
总结
至此,blob算法也就讲的差多不了,当然,这里面还有很多可以优化的空间。比如现在想到的是:链表的删除不用isDelete标志,而是真正的删除。blob特征计算(宽度,高度,面积,能量等)。游程结果按照16个一组拷贝到一个地方,使用AVX2等快速计算blob特征。游程提取步骤的并行加速。仿照FPGA的实现,对游程连接部分进行并行优化,等等。要优化的地方还很多,可惜假期要结束了,以后如果在工作中遇到了blob算法的优化问题,可能会接着往下做吧。
这次blob算法的调研与实现,更加深刻体会到了算法工程师查资料2天,编码2小时,调试2天的工作性质。作为算法工程师,一定不能闭门造车,要多方查阅资料,站在巨人的肩膀上。方案考虑周全,开始编码,真正编码时间肯定很短。接下来就是各种测试输入的异常调试,最后就是为了满足工程应用的各种简化与优化。