如何计算联通域重心

计算联通域重心

算法简介

     这是一幅阿米巴的图片，在假设质量平均的情况下也就是不考虑联通域加权的情况下我们可以使用“平衡法”来计算出它的重心。

     这是使用“平衡法”得到的计算结果，左边为原始图像，中间是经过二值化后再腐蚀10次的图像，最右边一幅使用蓝色圆圈标出了重心的坐标。

     为了对比方方便我们直接使用一个黑心圆继续测试程序，正如期望的那样，目测重心在圆心“附近”。

代码部分

     以下是代码部分，如果你已经配置好了OPENCV，应该可以跑起来。

     我只是用OPENCV来载入位图图像，图像的重心估算任务由函数 calCOG(...) 来完成。

// test.cpp 文件

// 计算联通域重心//#include "stdafx.h"#include <opencv.hpp>#include <iostream>#include "calCOG.hpp"using namespace std;using namespace cv;int _tmain(int argc, _TCHAR* argv[]){IplImage *src = cvLoadImage("data//circle.bmp");IplImage *result = cvCloneImage(src);cvNamedWindow("阿米巴原虫");cvShowImage("阿米巴原虫", src);// 对图像进行二值化cvThreshold(src, src, 200, 255, CV_THRESH_BINARY);// 腐蚀图像IplImage *eroderesult = cvCloneImage(src);cvErode(src, eroderesult, 0, 10);cvNamedWindow("腐蚀后的阿米巴");cvShowImage("腐蚀后的阿米巴",eroderesult);// 计算阿米巴的重心CvPoint cog;CvPoint anchor = cvPoint(250, 200);bool ok = calAmoebaCOG(eroderesult, anchor, cog);if( ok ){cout<<"阿米巴的重心坐标："<<"<"<<cog.x<<","<<cog.y<<">"<<endl;cvCircle(result, cog, 4, cvScalar(255,0,0, 2), 2);cvNamedWindow("阿米巴的重心");cvShowImage("阿米巴的重心",result);}else{cout<<"数据有误，计算失败！！"<<endl;}cvWaitKey();cvReleaseImage(&src);cvReleaseImage(&result);cvReleaseImage(&eroderesult);cvDestroyAllWindows();return 0;}

// calCOG.hpp 文件

#pragma once#include <iostream>#include <vector>#include <opencv.hpp>using namespace cv;using namespace std;// 如果越界，表达式值为false,否则为true#ifndef CHECKRANGE#define CHECKRANGE(value,min,max) (((value)>=(min))&&((value)<=(max)))#endifclass rowORcol{public:rowORcol();int id;// 行号或者列号int weight;// 当前行/列的重量int absdiff;// 被当前行/列分割的联通域两部分之差的绝对值};rowORcol::rowORcol():weight(0){}// 求两double之差绝对值double absdiff_double(const double &a,const double &b){double result;(a >= b)? (result = (a-b)):(result = (b-a));return result;}// 递归遍历联通域中所有像素（上下左右方向）// @copy:图像副本// @anchor_value:铆点像素第一通道值// @current:当前被遍历像素坐标// @rowdata:一个容器指针存放行统计数据// @coldata一个容器指针存放列统计数据void visitNeibor(IplImage *copy,const char anchor_value, CvPoint current, vector<rowORcol> *rowdata, vector<rowORcol> *coldata){char *current_ptr = copy->imageData + copy->widthStep*current.y/sizeof(char) + copy->nChannels*current.x;char current_value = *current_ptr;if( current_value == anchor_value ){if( current_value == 0 ){for(int x=0;x < copy->nChannels;++x){current_ptr[x] = 64;}}else{for(int x=0;x < copy->nChannels;++x){current_ptr[x] -= 1;}}bool findit = false;// 填充 rowdatafor(vector<rowORcol>::size_type x=0;x < rowdata->size();++x){if( current.y == rowdata->at(x).id ){findit = true;++(rowdata->at(x).weight);break;}}if( !findit ){rowORcol newrow;newrow.id = current.y;++newrow.weight;rowdata->push_back(newrow);}// 填充 coldatafindit = false;for(vector<rowORcol>::size_type x=0;x < coldata->size();++x){if( current.x == coldata->at(x).id ){findit = true;++(coldata->at(x).weight);break;}}if( !findit ){rowORcol newcol;newcol.id = current.x;++newcol.weight;coldata->push_back(newcol);}// 访问左邻居if( (current.x- 1) >= 0 ){CvPoint left = cvPoint(current.x-1,current.y);visitNeibor(copy, anchor_value, left, rowdata, coldata);}// 访问右邻居if( (current.y + 1) <= (copy->width -1) ){CvPoint right = cvPoint(current.x+1,current.y);visitNeibor(copy, anchor_value, right, rowdata, coldata);}// 访问上邻居if( (current.y - 1) >= 0 ){CvPoint up = cvPoint(current.x, current.y - 1);visitNeibor(copy, anchor_value, up, rowdata, coldata);}// 访问下邻居if( (current.y + 1) <+ (copy->height -1) ){CvPoint down = cvPoint(current.x, current.y + 1);visitNeibor(copy, anchor_value, down, rowdata, coldata);}}}// 计算铆点所属联通域的重心(COG,center_of_gravity)// @img:一幅二值图像指针// @anchor::铆点坐标，必须位于联通域内部// @cog:重心坐标// return value:成功返回true,否则返回falsebool calAmoebaCOG(const IplImage *img,const CvPoint anchor,CvPoint &cog){bool state = true;if( (!img)||!(CHECKRANGE(anchor.x,0,img->width))||!(CHECKRANGE(anchor.y,0,img->height))){state = false;}else{/* 算法简介 * 使用 calCenter(...) 得到联通域几何中心(anchor)  * if(weight_left(anchor.X) > (weight_right(anchor.x))) --anchor.x; * else ++anchor.x; * 达到平衡之后cog.x = anchor.x; * (计算重心Y坐标) * if(weight_up(anchor.y) > weight_down(anchor.y)) --anchor.y; * else --cnchor_y; * 达到平衡后cog.y = anchor.y; */IplImage *copy = cvCloneImage(img);vector<rowORcol> rowdata;vector<rowORcol> coldata;char anchor_value = *(img->imageData + img->widthStep*anchor.y/sizeof(char) + img->nChannels*anchor.x);CvPoint current = cvPoint(anchor.x, anchor.y);visitNeibor(copy, anchor_value, current, &rowdata, &coldata);                                ////////////////////////////////////////////////////////// 计算 weight.y// 将 rowdata 中元素按 id 进行升序排列vector<rowORcol>::size_type i,j;for(i=0;i < (rowdata.size() - 1);++i){rowORcol temp;for(j=0;j < (rowdata.size() - 1 - i);++j){if( rowdata.at(j).id > rowdata.at(j+1).id ){temp.id= rowdata.at(j+1).id;temp.weight= rowdata.at(j+1).weight;rowdata.at(j+1).id= rowdata.at(j).id;rowdata.at(j+1).weight= rowdata.at(j).weight;rowdata.at(j).id= temp.id;rowdata.at(j).weight= temp.weight;}}}// 计算 absdiffint rowsize = (int)rowdata.size();double *leftweight = new double[rowsize];double *rightweight = new double[rowsize];for(int x=0;x < rowsize;++x){if( x == 0 ){leftweight[0] = rowdata.at(0).weight;rightweight[rowsize-1] = rowdata.at(rowsize-1).weight;}else{leftweight[x] = leftweight[x-1] + (double)(rowdata.at(x).weight);rightweight[rowsize-1-x] = rightweight[rowsize-x] + (double)(rowdata.at(rowsize-1-x).weight);  }}for(vector<rowORcol>::size_type x=0;x < rowdata.size();++x){rowdata.at(x).absdiff = (int)absdiff_double(leftweight[x],rightweight[x]);}// 将 rowdata 中元素按 absdiff 进行升序排列for(i=0;i < (rowdata.size() - 1);++i){rowORcol temp;for(j=0;j < (rowdata.size() - 1 - i);++j){if( rowdata.at(j).absdiff > rowdata.at(j+1).absdiff ){temp.id= rowdata.at(j+1).id;temp.weight= rowdata.at(j+1).weight;temp.absdiff= rowdata.at(j+1).absdiff;rowdata.at(j+1).id= rowdata.at(j).id;rowdata.at(j+1).weight= rowdata.at(j).weight;rowdata.at(j+1).absdiff= rowdata.at(j).absdiff;rowdata.at(j).id= temp.id;rowdata.at(j).weight= temp.weight;rowdata.at(j).absdiff= temp.absdiff;}}}cog.y = rowdata.at(0).id;  // //////////////////////////////////////////////////////////////// 计算 weight.x// 计算每列的 absdiff// 升序排列// 首元素的id即为所求 weight.y// coldata 按 id(元素的列数)将元素升序排列for(i=0;i < (coldata.size() - 1);++i){rowORcol temp;for(j=0;j < (coldata.size() - 1 - i);++j){if( coldata.at(j).id > coldata.at(j+1).id ){temp.id= coldata.at(j+1).id;temp.weight= coldata.at(j+1).weight;coldata.at(j+1).id= coldata.at(j).id;coldata.at(j+1).weight= coldata.at(j).weight;coldata.at(j).id= temp.id;coldata.at(j).weight= temp.weight;}}}// 计算 absdiffint colsize = coldata.size();double *upweight = new double[coldata.size()];double *downweight = new double[coldata.size()];for(int x=0;x < colsize;++x){if( x == 0 ){upweight[0] = coldata.at(0).weight;downweight[colsize-1] = coldata.at(colsize-1).weight;}else{upweight[x] = upweight[x-1] + (double)(coldata.at(x).weight);downweight[colsize-1-x] = downweight[colsize-x] + (double)(coldata.at(colsize-1-x).weight);  }}for(vector<rowORcol>::size_type x=0;x < coldata.size();++x){coldata.at(x).absdiff = (int)absdiff_double(upweight[x],downweight[x]);}// 将 coldata 中元素按 absdiff 进行升序排列for(i=0;i < (coldata.size() - 1);++i){rowORcol temp;for(j=0;j < (coldata.size() - 1 - i);++j){if( coldata.at(j).absdiff > coldata.at(j+1).absdiff ){temp.id= coldata.at(j+1).id;temp.weight= coldata.at(j+1).weight;temp.absdiff= coldata.at(j+1).absdiff;coldata.at(j+1).id= coldata.at(j).id;coldata.at(j+1).weight= coldata.at(j).weight;coldata.at(j+1).absdiff= coldata.at(j).absdiff;coldata.at(j).id= temp.id;coldata.at(j).weight= temp.weight;coldata.at(j).absdiff= temp.absdiff;}}}cog.x = coldata.at(0).id;  // cvReleaseImage(©);delete []leftweight;delete []rightweight;delete []upweight;delete []downweight;}return state;}

     上边提供的“平衡法”并没有考虑到加权的问题，但是这是一种基本的方法，进行适当修改可以适合更多场合。

thx,that's all!

[e-mail]