图像处理找圆

2013-12-10 23:07:36 jia20003 阅读数 45444
  • 塑料瓶盖的度检测

    通过本课程学习,初学者可以熟练使用opencv4 API完成各种常见的图像分割、识别任务。 中高级学习者也一定能有新的体会和收获。 本课程所讲的例子代码来自于实际开发项目,有较高的实用性。

    18人学习 刘山
    免费试看

图像处理之霍夫变换圆检测算法

之前写过一篇文章讲述霍夫变换原理与利用霍夫变换检测直线, 结果发现访问量还是蛮

多,有点超出我的意料,很多人都留言说代码写得不好,没有注释,结构也不是很清晰,所以

我萌发了再写一篇,介绍霍夫变换圆检测算法,同时也尽量的加上详细的注释,介绍代码

结构.让更多的人能够读懂与理解.

一:霍夫变换检测圆的数学原理


根据极坐标,圆上任意一点的坐标可以表示为如上形式, 所以对于任意一个圆, 假设

中心像素点p(x0, y0)像素点已知, 圆半径已知,则旋转360由极坐标方程可以得到每

个点上得坐标同样,如果只是知道图像上像素点, 圆半径,旋转360°则中心点处的坐

标值必定最强.这正是霍夫变换检测圆的数学原理.

二:算法流程

该算法大致可以分为以下几个步骤


三:运行效果

图像从空间坐标变换到极坐标效果, 最亮一点为圆心.


图像从极坐标变换回到空间坐标,检测结果显示:


四:关键代码解析

个人觉得这次注释已经是非常的详细啦,而且我写的还是中文注释

	/**
	 * 霍夫变换处理 - 检测半径大小符合的圆的个数
	 * 1. 将图像像素从2D空间坐标转换到极坐标空间
	 * 2. 在极坐标空间中归一化各个点强度,使之在0〜255之间
	 * 3. 根据极坐标的R值与输入参数(圆的半径)相等,寻找2D空间的像素点
	 * 4. 对找出的空间像素点赋予结果颜色(红色)
	 * 5. 返回结果2D空间像素集合
	 * @return int []
	 */
	public int[] process() {

		// 对于圆的极坐标变换来说,我们需要360度的空间梯度叠加值
		acc = new int[width * height];
		for (int y = 0; y < height; y++) {
			for (int x = 0; x < width; x++) {
				acc[y * width + x] = 0;
			}
		}
		int x0, y0;
		double t;
		for (int x = 0; x < width; x++) {
			for (int y = 0; y < height; y++) {

				if ((input[y * width + x] & 0xff) == 255) {

					for (int theta = 0; theta < 360; theta++) {
						t = (theta * 3.14159265) / 180; // 角度值0 ~ 2*PI
						x0 = (int) Math.round(x - r * Math.cos(t));
						y0 = (int) Math.round(y - r * Math.sin(t));
						if (x0 < width && x0 > 0 && y0 < height && y0 > 0) {
							acc[x0 + (y0 * width)] += 1;
						}
					}
				}
			}
		}

		// now normalise to 255 and put in format for a pixel array
		int max = 0;

		// Find max acc value
		for (int x = 0; x < width; x++) {
			for (int y = 0; y < height; y++) {

				if (acc[x + (y * width)] > max) {
					max = acc[x + (y * width)];
				}
			}
		}

		// 根据最大值,实现极坐标空间的灰度值归一化处理
		int value;
		for (int x = 0; x < width; x++) {
			for (int y = 0; y < height; y++) {
				value = (int) (((double) acc[x + (y * width)] / (double) max) * 255.0);
				acc[x + (y * width)] = 0xff000000 | (value << 16 | value << 8 | value);
			}
		}
		
		// 绘制发现的圆
		findMaxima();
		System.out.println("done");
		return output;
	}
完整的算法源代码, 已经全部的加上注释

package com.gloomyfish.image.transform.hough;
/***
 * 
 * 传入的图像为二值图像,背景为黑色,目标前景颜色为为白色
 * @author gloomyfish
 * 
 */
public class CircleHough {

	private int[] input;
	private int[] output;
	private int width;
	private int height;
	private int[] acc;
	private int accSize = 1;
	private int[] results;
	private int r; // 圆周的半径大小

	public CircleHough() {
		System.out.println("Hough Circle Detection...");
	}

	public void init(int[] inputIn, int widthIn, int heightIn, int radius) {
		r = radius;
		width = widthIn;
		height = heightIn;
		input = new int[width * height];
		output = new int[width * height];
		input = inputIn;
		for (int y = 0; y < height; y++) {
			for (int x = 0; x < width; x++) {
				output[x + (width * y)] = 0xff000000; //默认图像背景颜色为黑色
			}
		}
	}

	public void setCircles(int circles) {
		accSize = circles; // 检测的个数
	}
	
	/**
	 * 霍夫变换处理 - 检测半径大小符合的圆的个数
	 * 1. 将图像像素从2D空间坐标转换到极坐标空间
	 * 2. 在极坐标空间中归一化各个点强度,使之在0〜255之间
	 * 3. 根据极坐标的R值与输入参数(圆的半径)相等,寻找2D空间的像素点
	 * 4. 对找出的空间像素点赋予结果颜色(红色)
	 * 5. 返回结果2D空间像素集合
	 * @return int []
	 */
	public int[] process() {

		// 对于圆的极坐标变换来说,我们需要360度的空间梯度叠加值
		acc = new int[width * height];
		for (int y = 0; y < height; y++) {
			for (int x = 0; x < width; x++) {
				acc[y * width + x] = 0;
			}
		}
		int x0, y0;
		double t;
		for (int x = 0; x < width; x++) {
			for (int y = 0; y < height; y++) {

				if ((input[y * width + x] & 0xff) == 255) {

					for (int theta = 0; theta < 360; theta++) {
						t = (theta * 3.14159265) / 180; // 角度值0 ~ 2*PI
						x0 = (int) Math.round(x - r * Math.cos(t));
						y0 = (int) Math.round(y - r * Math.sin(t));
						if (x0 < width && x0 > 0 && y0 < height && y0 > 0) {
							acc[x0 + (y0 * width)] += 1;
						}
					}
				}
			}
		}

		// now normalise to 255 and put in format for a pixel array
		int max = 0;

		// Find max acc value
		for (int x = 0; x < width; x++) {
			for (int y = 0; y < height; y++) {

				if (acc[x + (y * width)] > max) {
					max = acc[x + (y * width)];
				}
			}
		}

		// 根据最大值,实现极坐标空间的灰度值归一化处理
		int value;
		for (int x = 0; x < width; x++) {
			for (int y = 0; y < height; y++) {
				value = (int) (((double) acc[x + (y * width)] / (double) max) * 255.0);
				acc[x + (y * width)] = 0xff000000 | (value << 16 | value << 8 | value);
			}
		}
		
		// 绘制发现的圆
		findMaxima();
		System.out.println("done");
		return output;
	}

	private int[] findMaxima() {
		results = new int[accSize * 3];
		int[] output = new int[width * height];
		
		// 获取最大的前accSize个值
		for (int x = 0; x < width; x++) {
			for (int y = 0; y < height; y++) {
				int value = (acc[x + (y * width)] & 0xff);

				// if its higher than lowest value add it and then sort
				if (value > results[(accSize - 1) * 3]) {

					// add to bottom of array
					results[(accSize - 1) * 3] = value; //像素值
					results[(accSize - 1) * 3 + 1] = x; // 坐标X
					results[(accSize - 1) * 3 + 2] = y; // 坐标Y

					// shift up until its in right place
					int i = (accSize - 2) * 3;
					while ((i >= 0) && (results[i + 3] > results[i])) {
						for (int j = 0; j < 3; j++) {
							int temp = results[i + j];
							results[i + j] = results[i + 3 + j];
							results[i + 3 + j] = temp;
						}
						i = i - 3;
						if (i < 0)
							break;
					}
				}
			}
		}

		// 根据找到的半径R,中心点像素坐标p(x, y),绘制圆在原图像上
		System.out.println("top " + accSize + " matches:");
		for (int i = accSize - 1; i >= 0; i--) {
			drawCircle(results[i * 3], results[i * 3 + 1], results[i * 3 + 2]);
		}
		return output;
	}

	private void setPixel(int value, int xPos, int yPos) {
		/// output[(yPos * width) + xPos] = 0xff000000 | (value << 16 | value << 8 | value);
		output[(yPos * width) + xPos] = 0xffff0000;
	}

	// draw circle at x y
	private void drawCircle(int pix, int xCenter, int yCenter) {
		pix = 250; // 颜色值,默认为白色

		int x, y, r2;
		int radius = r;
		r2 = r * r;
		
		// 绘制圆的上下左右四个点
		setPixel(pix, xCenter, yCenter + radius);
		setPixel(pix, xCenter, yCenter - radius);
		setPixel(pix, xCenter + radius, yCenter);
		setPixel(pix, xCenter - radius, yCenter);

		y = radius;
		x = 1;
		y = (int) (Math.sqrt(r2 - 1) + 0.5);
		
		// 边缘填充算法, 其实可以直接对循环所有像素,计算到做中心点距离来做
		// 这个方法是别人写的,发现超赞,超好!
		while (x < y) {
			setPixel(pix, xCenter + x, yCenter + y);
			setPixel(pix, xCenter + x, yCenter - y);
			setPixel(pix, xCenter - x, yCenter + y);
			setPixel(pix, xCenter - x, yCenter - y);
			setPixel(pix, xCenter + y, yCenter + x);
			setPixel(pix, xCenter + y, yCenter - x);
			setPixel(pix, xCenter - y, yCenter + x);
			setPixel(pix, xCenter - y, yCenter - x);
			x += 1;
			y = (int) (Math.sqrt(r2 - x * x) + 0.5);
		}
		if (x == y) {
			setPixel(pix, xCenter + x, yCenter + y);
			setPixel(pix, xCenter + x, yCenter - y);
			setPixel(pix, xCenter - x, yCenter + y);
			setPixel(pix, xCenter - x, yCenter - y);
		}
	}

	public int[] getAcc() {
		return acc;
	}

}
测试的UI类:

package com.gloomyfish.image.transform.hough;

import java.awt.BorderLayout;
import java.awt.Color;
import java.awt.Dimension;
import java.awt.FlowLayout;
import java.awt.Graphics;
import java.awt.Graphics2D;
import java.awt.GridLayout;
import java.awt.event.ActionEvent;
import java.awt.event.ActionListener;
import java.awt.image.BufferedImage;
import java.io.File;

import javax.imageio.ImageIO;
import javax.swing.BorderFactory;
import javax.swing.JButton;
import javax.swing.JFrame;
import javax.swing.JPanel;
import javax.swing.JSlider;
import javax.swing.event.ChangeEvent;
import javax.swing.event.ChangeListener;

public class HoughUI extends JFrame implements ActionListener, ChangeListener {
	/**
	 * 
	 */
	public static final String CMD_LINE = "Line Detection";
	public static final String CMD_CIRCLE = "Circle Detection";
	private static final long serialVersionUID = 1L;
	private BufferedImage sourceImage;
// 	private BufferedImage houghImage;
	private BufferedImage resultImage;
	private JButton lineBtn;
	private JButton circleBtn;
	private JSlider radiusSlider;
	private JSlider numberSlider;
	public HoughUI(String imagePath)
	{
		super("GloomyFish-Image Process Demo");
		try{
			File file = new File(imagePath);
			sourceImage = ImageIO.read(file);
		} catch(Exception e){
			e.printStackTrace();
		}
		initComponent();
	}
	
	private void initComponent() {
		int RADIUS_MIN = 1;
		int RADIUS_INIT = 1;
		int RADIUS_MAX = 51;
		lineBtn = new JButton(CMD_LINE);
		circleBtn = new JButton(CMD_CIRCLE);
		radiusSlider = new JSlider(JSlider.HORIZONTAL, RADIUS_MIN, RADIUS_MAX, RADIUS_INIT);
		radiusSlider.setMajorTickSpacing(10);
		radiusSlider.setMinorTickSpacing(1);
		radiusSlider.setPaintTicks(true);
		radiusSlider.setPaintLabels(true);
		numberSlider = new JSlider(JSlider.HORIZONTAL, RADIUS_MIN, RADIUS_MAX, RADIUS_INIT);
		numberSlider.setMajorTickSpacing(10);
		numberSlider.setMinorTickSpacing(1);
		numberSlider.setPaintTicks(true);
		numberSlider.setPaintLabels(true);
		
		JPanel sliderPanel = new JPanel();
		sliderPanel.setLayout(new GridLayout(1, 2));
		sliderPanel.setBorder(BorderFactory.createTitledBorder("Settings:"));
		sliderPanel.add(radiusSlider);
		sliderPanel.add(numberSlider);
		JPanel btnPanel = new JPanel();
		btnPanel.setLayout(new FlowLayout(FlowLayout.RIGHT));
		btnPanel.add(lineBtn);
		btnPanel.add(circleBtn);
		
		
		JPanel imagePanel = new JPanel(){

			private static final long serialVersionUID = 1L;

			protected void paintComponent(Graphics g) {
				if(sourceImage != null)
				{
					Graphics2D g2 = (Graphics2D) g;
					g2.drawImage(sourceImage, 10, 10, sourceImage.getWidth(), sourceImage.getHeight(),null);
					g2.setPaint(Color.BLUE);
					g2.drawString("原图", 10, sourceImage.getHeight() + 30);
					if(resultImage != null)
					{
						g2.drawImage(resultImage, resultImage.getWidth() + 20, 10, resultImage.getWidth(), resultImage.getHeight(), null);
						g2.drawString("最终结果,红色是检测结果", resultImage.getWidth() + 40, sourceImage.getHeight() + 30);
					}
				}
			}
			
		};
		this.getContentPane().setLayout(new BorderLayout());
		this.getContentPane().add(sliderPanel, BorderLayout.NORTH);
		this.getContentPane().add(btnPanel, BorderLayout.SOUTH);
		this.getContentPane().add(imagePanel, BorderLayout.CENTER);
		
		// setup listener
		this.lineBtn.addActionListener(this);
		this.circleBtn.addActionListener(this);
		this.numberSlider.addChangeListener(this);
		this.radiusSlider.addChangeListener(this);
	}
	
	public static void main(String[] args)
	{
		String filePath = System.getProperty ("user.home") + "/Desktop/" + "zhigang/hough-test.png";
		HoughUI frame = new HoughUI(filePath);
		// HoughUI frame = new HoughUI("D:\\image-test\\lines.png");
		frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
		frame.setPreferredSize(new Dimension(800, 600));
		frame.pack();
		frame.setVisible(true);
	}

	@Override
	public void actionPerformed(ActionEvent e) {
		if(e.getActionCommand().equals(CMD_LINE))
		{
			HoughFilter filter = new HoughFilter(HoughFilter.LINE_TYPE);
			resultImage = filter.filter(sourceImage, null);
			this.repaint();
		}
		else if(e.getActionCommand().equals(CMD_CIRCLE))
		{
			HoughFilter filter = new HoughFilter(HoughFilter.CIRCLE_TYPE);
			resultImage = filter.filter(sourceImage, null);
			// resultImage = filter.getHoughSpaceImage(sourceImage, null);
			this.repaint();
		}
		
	}

	@Override
	public void stateChanged(ChangeEvent e) {
		// TODO Auto-generated method stub
		
	}
}
五:霍夫变换检测圆与直线的图像预处理

使用霍夫变换检测圆与直线时候,一定要对图像进行预处理,灰度化以后,提取

图像的边缘使用非最大信号压制得到一个像素宽的边缘, 这个步骤对霍夫变

换非常重要.否则可能导致霍夫变换检测的严重失真.

第一次用Mac发博文,编辑不好请见谅!

2016-04-26 16:46:27 jsgaobiao 阅读数 24868
  • 塑料瓶盖的度检测

    通过本课程学习,初学者可以熟练使用opencv4 API完成各种常见的图像分割、识别任务。 中高级学习者也一定能有新的体会和收获。 本课程所讲的例子代码来自于实际开发项目,有较高的实用性。

    18人学习 刘山
    免费试看

代码下载:http://download.csdn.net/detail/jsgaobiao/9503229

 

Ø  【作业要求】


Write your own imfindcircles() to simulatematlab function imfindcircles(). The attached images are for testing.

Submit your code, result and report.


Ø  【文件说明】


main.m:

读取3张测试图片并且进行转灰度图、用Sobel算子边缘检测、二值化、调用findcircle函数并将找到的圆绘制在原图上。

findcircle.m:

利用霍夫变换的方法找到图片中的圆,并返回圆的坐标和半径。


Ø  【作业思路】


本次作业找圆主要是使用了霍夫变换的方法。

先介绍一下霍夫变换的思想:


霍夫变换是图像处理中用于识别几何图形的一种常用方法。最简单的应用是检测黑白图像上的一条直线。我们用y=kx+b表示一条直线,其中k和b是参数,分别是斜率和截距。如果将k和b视为变量,考虑参数平面k-b,那么x-y平面上的一个点对应到k-b参数平面上就是一条直线。

假设在x-y平面上有若干待拟合的散点,那么他们可以对应到参数平面上的若干条直线。如果x-y平面上的散点都在一条直线上,那么他们对应在参数平面上的若干条直线必然交与(k0,b0)这一点,其中k0,b0是x-y平面上那条直线的参数。如果参数平面上的直线产生了多个交点,那么我们就选取经过直线最多的交点。可以理解为若干直线对于坐标点(k,b)的投票,票数最多的参数(k,b)我们认为就是最好的一组参数。


处理圆的情况更容易理解,因为我们可以取和图像平面一样的参数平面,以图像上的每一个非零点(因为该算法处理的是二值化后的图片,即对原图进行边缘提取后的结果)为圆心,以已知的半径在参数平面上画圆(对圆覆盖的坐标点进行投票),最后找出参数平面上的峰值,就对应于原图中的圆心。

但是这个做法存在一个问题——我们并不知道圆的半径。一个解决方案是把参数平面拓展为三维的空间x-y-r,对应圆心坐标和半径。这样做需要增加相当大的时间、空间复杂度。为了提高算法的效率,我将圆的半径范围作为参数输入findcircle函数中,对于不同的图像采用不同的参数进行处理,可以有效的减少运行的时间。

 

另外,作为参数输入findcircle函数的还有二值化后的图像BW,枚举半径的范围[minR,maxR]和步长stepR和在参数平面画圆投票的弧度步长stepAngle。

最终的检测效果如下图所示:

2019-03-16 14:32:13 weixin_44540503 阅读数 1750
  • 塑料瓶盖的度检测

    通过本课程学习,初学者可以熟练使用opencv4 API完成各种常见的图像分割、识别任务。 中高级学习者也一定能有新的体会和收获。 本课程所讲的例子代码来自于实际开发项目,有较高的实用性。

    18人学习 刘山
    免费试看

利用opencv进行图像处理,提取椭圆圆心处理

写这个是因为项目正好在做这个,所以简单写写提取椭圆圆心坐标的代码,用的软件是VS。
首先介绍一下步骤,直接从图像处理开始
1,二值化处理(threhold())
2,高斯滤波(GaussianBlur())
3,轮廓提取(canny算子)
4,寻找闭合轮廓(findContours())

上面介绍的是提取椭圆轮廓过程,下面介绍的是过滤干扰图像的过程,分三步,不过,这三部也不都是必须的,看个人需求:
1,像素点数量过滤(size)
2,面积过滤
3,长宽比过滤
代码如下,比较简单,就不过多介绍了,欢迎留言:

#include<iostream>
#include"opencv.hpp"

using namespace std;
using namespace cv;


void drawCross(Mat &img, Point2f point, Scalar color, int size, int thickness /*= 1*/)
{
	//绘制横线  
	line(img, cvPoint(point.x - size / 2, point.y), cvPoint(point.x + size / 2, point.y), color, thickness, 8, 0);
	//绘制竖线  	
	line(img, cvPoint(point.x, point.y - size / 2), cvPoint(point.x, point.y + size / 2), color, thickness, 8, 0);
	return;
}
 

//提取单幅图像的特征点
void FeaturePoint(Mat &img)
{
	Point2f center; //定义变量
	vector<Point2f> ellipsecenterleft;
	Mat edges;
	threshold(img, img, 80, 255, CV_THRESH_BINARY_INV);
	//imshow("threshold", img);
	////此函数等待按键,按键盘任意键就返回
	//waitKey(0);
	GaussianBlur(img, edges, Size(5, 5), 0, 0);
	//imshow("GaussianBlur", edges);
	////此函数等待按键,按键盘任意键就返回
	//waitKey(0);
	Canny(edges, edges, 40, 120, 3);
	//imshow("Canny", edges);
	////此函数等待按键,按键盘任意键就返回
	//waitKey(0);
	vector<vector<Point> > contours;// 创建容器,存储轮廓
	vector<Vec4i> hierarchy;// 寻找轮廓所需参数


	findContours(edges, contours, hierarchy, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE);
	//cout << "conours size in "<<serialNumber <<" :" << contours.size() << endl;

	//if (contours.size() == 1)
	for (int i = 0; i < contours.size();i++) 
	{
		RotatedRect m_ellipsetemp;  // fitEllipse返回值的数据类型
    	if (contours[i].size() <= 200){
			continue;
		}
		if (contourArea(contours[i]) < 100 && contourArea(contours[i]) > 1000)
		{
			continue;
		}
		m_ellipsetemp = fitEllipse(contours[i]);  //找到的第一个轮廓,放置到m_ellipsetemp


		ellipse(img, m_ellipsetemp, cv::Scalar(255,0,0));   //在图像中绘制椭圆

		if (m_ellipsetemp.size.width / m_ellipsetemp.size.height <0.3)
		{
			continue;
		}
		center = m_ellipsetemp.center;//读取椭圆中心

		drawCross(img, center, Scalar(255,0,0), 30, 2);
		cout << center.x << ends << center.y << endl;
	}
	
	/*imshow("image", img);
	waitKey(0);*/
	//return center;//返回椭圆中心坐标
	
}



int main(int argc, char* argv[])
{
	const char* imagename = "pixmap1.png";

	//从文件中读入图像
	Mat img = imread(imagename, 1);

	//如果读入图像失败
	if (img.empty())
	{
		fprintf(stderr, "Can not load image %s\n", imagename);
		return -1;
	}
	FeaturePoint(img);
	//显示图像
	imshow("image", img);

	//此函数等待按键,按键盘任意键就返回
	waitKey();
	return 0;
}
2019-05-16 10:11:18 cheng_xing_ 阅读数 1920
  • 塑料瓶盖的度检测

    通过本课程学习,初学者可以熟练使用opencv4 API完成各种常见的图像分割、识别任务。 中高级学习者也一定能有新的体会和收获。 本课程所讲的例子代码来自于实际开发项目,有较高的实用性。

    18人学习 刘山
    免费试看

经典椭圆检测方法

椭圆检测算法经过多年的研究发展,已经基本形成一个较完整的体系。它们大致可以分为三类即投票(聚类)、最优化、基于弧段的方法。

投票(聚类)方法

椭圆因为有中心位置坐标、长短轴长度、倾斜角五个参数,标准霍夫变换有较强的鲁棒性,但对内存要求高,运算效率低,不太现实。霍夫变换类算法以霍夫变换为算法基础,经过不同国家研究人员多年的不懈努力研究,如今已衍生出很多改进算法,它们各有优劣。随机霍夫变换算法相对标准霍夫变换计算速度有较大提升,但检测相互遮挡的椭圆时准确度低。

随机hough变换椭圆检测算法

随机椭圆检测结合使用了了最小二乘法和Hough变换算法。第一步预处理,获得较理想的边缘图。第二步随机选取三个点,取这三点为中心相同大小的邻域中所有点,用最小二乘法把它们拟合成一个椭圆。如图2-3所示。第三步从边缘点中再随机选取第四个点,判断此点是否在拟合出的圆上。若是,则认为该椭圆是真实椭圆的可能性较大,接着收集证据,验证该椭圆的真实性。

图2-3 随机选点示意图
算法具体过程如下(从第二步开始):
1.把边缘检测得到的点收进集合V中,失败计数器f初始值设为0。设定5个阈值,分别是能容忍的失败次数最大值Tf,检测进行时对V中边缘点数量的要求阈值Tem,随机选取的三点之间两两距离最小值Ta,随机选取的第四点到可能椭圆边界距离的最大值Td,以及椭圆残缺比率阈值Tr。
2.np表示集合V中剩余的点的数量,当np小于Tem时或当失败次数f大于Tf时停止检测,算法终止;否则从V中随机取四点,并从V中删除这四点。
3.若用来求解椭圆参数的三个点两两之间距离都大于Ta,拟合出椭圆,计算第四个点到该椭圆边界的距离,若距离小于Td,执行第4步;若不满足两者之一,将这四个点返回到V中,失败次数加一,回到第2步执行。
4.设E为第3步拟合出来的椭圆,初始化满足阈值Td的点的个数num。遍历V中点,计算并判断它们到椭圆E的边界的距离是否小于Td,若是则num=num+1,并将该点从V中除去,直到遍历完成。
5.若num>=Tr*K,其中K为椭圆E的周长,那么跳转到第6步;否则认为椭圆E不是真实的椭圆,将第4步和第2步中删除的num+4个点返回V中,并跳转到第2步。
6.认为椭圆E是一个真实存在的椭圆,f置0,并跳转到第2步。
随机hough变换的优缺点如下:
第一,由于该算法是基于最小二乘法,所以一方面检测结果往往比真实椭圆小而且对噪声敏感,但是另一方面当预处理效果较好时检测精度很高。
第二,由于该算法基于随机采样,所以一方面可能会有所选点距离较近的情况造成拟合出的椭圆偏差较大,但是另一方面因为随机采样的灵活性检测速度提升了。
第三,一方面当参数选取的较好时检测又快有准确;另一方面,由于该算法严格由参数Ta,Td,Tr控制而且这些参数不易取到合适值,所以会出现不合适的参数不仅增加计算量,而且增加误检机会的情况。

最优化方法

最优化类方法优点在其精度上,缺点是其一次只能处理一个图形,即此前要对图像信息进行分类分离。AndrewFitzgibb等人提出了直接最小二乘法椭圆拟合算法。该方法能保证拟合出来的一定是椭圆。但该方法受到孤立点和噪声点的影响。目前最优化算法多与其他算法一起结合使用。

基于弧段的方法

基于边界聚类的椭圆检测方法结合使用了基于弧段的方法和最小二乘法。从边界图提取圆弧,再经过过滤、聚类,最终用最小二乘法拟合出椭圆。该方法能有效应对多个椭圆、椭圆相互遮挡和椭圆部分缺损等复杂情况,因而引起了广泛的注意。

边界聚类算法流程

边界聚类算法属于从下往上结构的算法。算法步骤主要分为三步,分别是预处理,边界聚类和直接最小二乘法拟合椭圆三个过程。流程图如下所示:
在这里插入图片描述

预处理

预处理第一步是进行灰度变换。灰度化常用的方法也就是依据亮度方程来实现的,即依据人眼对不同颜色的敏感度不同,对RGB分量以不同系数的加权平均。
在这里插入图片描述
第二是降噪。去噪手段对应于噪声的两种分类主要有两种。噪声功率谱符合高斯函数时用可以用高斯平滑模板平滑。由于脉冲干扰等产生的噪声(即椒盐噪声)可以采用中值滤波去除。
第三步是边界检测。通常图像中边界点都是图像中亮度梯度比较大的点,这些点包含了我们检测要用到的图像特征信息。边界检测最常用的方法是Canny算法。该算法主要分四个部分。一、降噪。方法是让原始图像和所用的高斯模板作卷积,模板在使用前指定了标准差。 二、寻找亮度梯度。Canny算法使用4个模板检测边缘的方向,,它们分别是水平、垂直、主对角线和副对角线方向;遍历圆图上的每个像素点,让原始图像中以该点为中心以该模板为窗口内的所有点与该模板作卷积,我们就从原图获得了各个点亮度梯度图和亮度梯度的方向。三、非极大值抑制。该过程目的是获得单像素的候选边缘,主要操作是将非零像素点所在的区域进行细化。具体过程如下:对于图3-2中一点P(x,y),计算P点梯度方向与其8-连通邻域点所组成的正方形的交点 (x1,y1)和(x2,y2)。如图2-5,交点坐标通过插值法得到。如果中间点的值大于这两个交点值,那么P点值不变,如若不然置零。四、滞后阈值操作。它需要设置两个阈值t1与t2。t1等于边界像素数除以总像素数,这些点称之为强边缘像素。t2等于t1除以2, t2和 t1之间的点称之为弱边缘像素。最后通过将8-连通的弱像素集成到强像素,再把它们连接起来,得到边界图。
在这里插入图片描述
第四步是二值化。二值化的效果几乎完全取决于分割阈值的选取。所以自动寻找最佳分割阈值的方法就显得十分关键。找到图片二值化的一个合适的分割阈值的一种方法是按图像的灰度特性,将图像分成背景和目标两部分,背景和目标之间的类间方差最大的分割意味着错分概率最小。MATLAB 的Graythresh函数就是使用该方法来获得一个自适应阈值作为二值化的分割依据。
第五步,二值化后,因为椭圆弧附近的非相关像素会严重影响检测结果,因此为了大幅度减少非相关像素,本文用形态学的腐蚀操作来得到细化的边界。

边界像素连接

采用Kovesi的边界连接算法,以8-邻域连通准则从上至下,从左至右扫描二值图像,将边界像素连接为有向边界列。然后采用边界列中像素数阈值条件去除像素数较小的集合。因为若边界列像素数少于阈值数,则很有可能是噪声或背景,应当删除。具体步骤如下:
1.以8-邻域连通准则从上至下,从左至右扫描二值图像,按连通区域对图像中的像素点聚类。
2.寻找每一个连通域中边界像素中所有的结束点和分叉点(分叉点是三条以上曲线的交点)并存储。
3.以这些结束点和分叉点为结束标志,让每一个连通域中的像素点集合分割为遇到结束点和连接点就断开的小集合。
4.删除这些集合中像素数小于某一阈值的部分。

线段列提取

因为图像光栅化难以获得准确的切线,而后续过程需要用到圆弧的切线,所以要进行线段拟合,即用多段折线代替原来的圆弧。具体步骤如下:
1.取边界像素连接成的第i条有向边界列,判断是否超过边界列总数目total,若不是进行步骤2,若是终止算法。
2.判断该边界列是否已经完成处理,若为否则进行第3步;若是则i=i+1,重新进行步骤2.
3.从其中第三个点开始,计算第一个点到这个点(记为点j)的连线方程,并依次判断第一个点和该点之间的所有点到该连线的距离,若所有距离均小于某一阈值,则j=j+1,重新进行步骤3,否则该有向边界列从该处断开,前面部分只保留第一个点和第j个点,前面j个点构成的连线用第一点和第j点之间的直线连线代替;后面部分仍然记为边界列i,
4.判断步骤3中断开的有向边界列后面部分像素数是否小于某阈值,若是则删除掉,否则不处理。最后跳转到步骤1。
完成这一个过程后一个连通域的的像素点构成的曲线就变成了其中部分像素点构成的一条折线。经过这个过程虽然像素信息损失了一部分,但是求取圆弧切线的精度从某种意义上说提高了,因为没有了光栅化效应。而且数据少了处理变得简单。再采用线段数阈值条件去除较短的线段列。若线段数数少于阈值数,则很有可能是噪声或背景或者进行拟合时误差过大,因此须删除。

线段列旋转方向统一

本文将所有线段列旋转方向统一为逆时针方向。
假设图3-4中的黑点为线段列中的点,箭头代表线段列的方向,P1(x1,y1),P2(x2,y2),P3(x3,y3)为线段列中连续的三个像素,像素都引入z坐标,且令其为0,则P1(x1,y1,0),P2(x2,y2,0),P3(x3,y3,0),空间向量
P1P2=(x2-x1,y2-y1,0)         (3-2)P2P3=(x3-x2,y3-y2,0)         (3-3)
向量积
P1P2×P2P3=|■(i ⃗&j ⃗&k ⃗@x2-x1&y2-y1&0@x3-x2&y3-y2 &0)|=(0,0,(x2-x1)(y3-y2)-(x3-x2)(y2-y1))         (3-4)
对一个线段列中除去首尾两个点的所有点像P2点一样计算并判别和存储,若小于0的次数最多,则认为线段列的旋转方向是顺时针,将线段列中的点逆序处理;若大于0的次数最多,则认为线段列的旋转方向是逆时针。
在这里插入图片描述

凹点和角点检测

在确定线段列的旋转方向为逆时针方向后,检测凹点和角点方法同前面统一线段列旋转方向类似,对一个线段列中除去首尾两个点的所有点计算P1P2×P2P3并判断向量积第三个分量的大小,若小于0,则P2为凹点。因为边界波动可能引入冗余凹点也即因边界检测误差可能错判一些正常点为凹点而进修分割会导致检测率下降,所以增加一个角度判断过程,即前面向量积为0并且P1P2和P2P3的夹角大于某阈值,才为凹点,这样选择合适的阈值可保证凹点检测的准确性。若向量积大于0且P1P2和P2P3的夹角大于另一阈值,则认为该点角度变化过大,是角点,线段列有很大可能性不是椭圆弧,而有可能是三角形、矩形等图形的边角,因此从该点分割该线段列。分割完成后,过滤掉含点数较少的线段列,即可除掉部分非椭圆弧。留下的线段列认为是椭圆弧,参加后续的聚类。如图3-5所示,左上角的P2很可能是凹点,右下角的P2很可能是角点。
在这里插入图片描述

圆弧聚类

圆弧聚类是将属于同一椭圆但是分开的两条或多条椭圆弧进行聚类。在进行聚类前,首先要判断弧段的完整度。一般用弧段的首尾端点P1,P3与中点P2构成的两向量P2P1,P2P3的夹角的大小来进行判断。夹角越小,一般该椭圆弧越完整,夹角越大,一般认为椭圆弧缺损越严重。设定一个阈值,当夹角小于该阈值时认为该弧段已经足够完整,仅仅靠该弧段上的点就可以较准确地拟合出真实存在的椭圆,因此该弧段不需要参与后面的聚类过程。如果希望该阈值自适应,在划分待聚类椭圆弧和不须聚类的弧(较完整弧)之前,先要确定该阈值。用直接最小二乘法拟合该弧所在的椭圆,若较圆,为了减小的误差,应使阈值夹角稍微大一些,如90度;若该弧所在的椭圆较扁,应使阈值夹角稍微小一些,如60度。接下来才根据该自适应阈值进行对圆弧判断。当大于阈值时认为该弧段上的点过少,不足以拟合出准确的椭圆,需要找到和该弧段属于同一椭圆的弧段然后用它们所有的点一起拟合出一个椭圆。经过此判断过程,椭圆弧就被分成两组。把需要参加聚类的椭圆弧按照含点数的数目由多到少进行排列,下面的过程都按照数目多的弧段优先的顺序进行。
对于待聚类的椭圆弧,先要定义其搜索区域,由于椭圆是封闭图形,所以整个椭圆可以确定是在其任何一部分弧和弧两端点的切线所在的射线包围起来的区域里面,属于该椭圆的其他弧以确定是在该弧对应的弦和弧两端点的切线所在的射线包围起来的区域里面。这就是我们搜索的区域。在图中a1的搜索区域也就是射线l1,l2,弦l3和图像边缘范围内的区域,在这个区域里面找弧,可以缩小搜寻的范围,提高效率。判断一条弧是否在待聚类椭圆弧的搜索区域里面我们只需取这条弧的首末端点j3,j4是否在搜索区域。方法是分别求过这两点同时平行于待聚类椭圆弧对应弦l3的直线和切线的交点,若交点分别有两个,交点都在射线上且这两个端点在对应两个交点之间则该弧段在搜索区域内。图中明显a2,a3,a4在a1的搜索区域内而a5不在。
在这里插入图片描述
待聚类椭圆弧找到待配对的圆弧后用两种约束条件判断它们到底是否属于同一椭圆。约束一是利用一个椭圆任意两段弧弧中点之间的距离大于一个弧中点到另一个弧首末端点连线的中点之间的距离,在图中即为
在这里插入图片描述
用来去除图中a2类型的椭圆弧。右图(b)不满足,左图(a)同时满足这两个关系,进入下一步,再用约束二进行判断。
在这里插入图片描述
约束二是点到拟合椭圆边界距离条件。我们只需要让两条线段列中的点一起参与椭圆拟合,按照下面公式计算所有这些点到拟合出椭圆边界的距离。设置一距离阈值,当d_i小于该阈值认为该点落在该椭圆上,否则该点不在这个椭圆上。统计d_i中小于某一阈值的点的数量,若大于某一比例(比例阈值),则认为这两条弧属于同一椭圆,否则不属于同一椭圆。判定后将属于同一椭圆的弧段聚类到一起。
在这里插入图片描述
其中
其中:x^'=(x_i-x_0)cos⁡〖θ+(y_i-y_0)sin⁡θ 〗,y^'=-(x_i-x_0)sin⁡〖θ+(y_i-y_0)cos⁡θ 〗。
在这里插入图片描述
如上图(a)中的参与拟合的点较多都落在拟合的椭圆上,所以有较大可能满足约束二条件;(b)中大多数拟合点点离拟合出的椭圆边界有一定距离,有较大可能不满足约束条件二。
如果希望这两个阈值改为自适应的,方法是先拟合出椭圆,判断椭圆大小。若椭圆较小,应适当降低限制,即增大距离阈值,减小比例阈值;若椭圆较大,应适当提高限制,即减小距离阈值,增大比例阈值。本算法中选取的是,若椭圆的短轴小于50则距离阈值为0.05,比例阈值设为0.7;否则前者取0.03,后者取0.8。

再配对

聚类后的弧和较完整弧或者两个较完整弧可能属于同一椭圆但是在前面的步骤它们只是被分开了并没有配对,所以有必要增加再匹配过程,增加检测准确度。匹配方法还是约束条件二的方法。因为该方法和原算法的去伪过程相似,所以经过该方法后无需再去伪。

直接最小二乘法椭圆拟合

前面的很多步骤都删除了像素点较少的集合,或者用少数点代替了边界列中的很多点,或者是分割后再删除点数较少的集合的,这些操作到椭圆拟合这一步实际上基本上去除了所有的背景和噪声,甚至包括一部分有用信息。所以即使对噪声和孤立点敏感的直接最小二乘法也可以用来拟合椭圆,而且因为该方法对椭圆缺损不敏感,所以非常适合。

实验效果

边界聚类算法检测结果

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

边界聚类算法和随机霍夫变换算法比较

用时比较在这里插入图片描述
检测结果比较
左图是选点情况;右图是拟合的椭圆和真实椭圆的差异
在这里插入图片描述
随机霍夫变换算法检测结果
在这里插入图片描述
边界聚类算法检测结果

2014-03-30 21:18:14 zhiyuanzhe007 阅读数 9368
  • 塑料瓶盖的度检测

    通过本课程学习,初学者可以熟练使用opencv4 API完成各种常见的图像分割、识别任务。 中高级学习者也一定能有新的体会和收获。 本课程所讲的例子代码来自于实际开发项目,有较高的实用性。

    18人学习 刘山
    免费试看

        图像处理中的椭圆检测用处还是挺多的,找到这里来的同学大多是想用椭圆检测来解决某些实际问题吧,所以我就不做介绍,直奔主题。我研究这块也有一段时间了,也查找了挺多资料,貌似通用的椭圆算法还没有,不排除我孤陋寡闻了。前辈提出的算法适用范围比较有限,这个比较有限是相对直线检测来说的。但直接用Hough变换来找椭圆几乎是不可能的事,在5维空间里做投票,想想都觉得可怕。于是有人想到用随机Hough变换。这是一种很合理的方法,我就是这么做的,不过这种方法有个不足之处,后面会讲到。这里先介绍这方法的流程。

        二次曲线的一般方程为:,其中(x,y)为图像的坐标空间,BCDEF是二次曲线的参数。当满足时二次曲线为椭圆。方程中需要求解的参数有5个,在随机Hough变换过程中要至少采集5个点,得到5个方程以求得这5个参数。若方程有解且满足约束条件,则将解加入参数空间进行累积。思路是比较简单的,下面边贴代码边解释(P.S. 代码仅供参考)。

void FindEllipse(TImage* OrgGrayImg)
{
int ImgHeight = OrgGrayImg->nHeight;
	int ImgWidth = OrgGrayImg->nWidth;
	unsigned char * Img = OrgGrayImg->pImage; // 输入图像确保是二值图像

srand((unsigned)time(NULL));

	int totalPt = 0;// 用于统计样本点的个数

	for (i = 0; i < ImgHeight - 0; i++)
	{
		unsigned char *imgdata = Img + i * ImgWidth;
		for (j = 0; j < ImgWidth - 0; j++)
		{
			if (!imgdata[j])
				totalPt ++;
		}
	}

	if (totalPt < 5)
		return;

	POINT * seq;
	seq = new POINT [totalPt];

	int count = 0;
	for (i = 0; i < ImgHeight; i++)
	{
		unsigned char *data = Img + i * ImgWidth;
		for (j = 0; j < ImgWidth; j++)
		{
			if (!data[j])
			{
				seq[count].x = j;
				seq[count].y = i;
				count ++;
			}
		}
	}

	double Para[5];	// 存放结果(5个参数A,B,C,D,E)的数组
	int Angle_V[360]={0};// 椭圆倾斜角参数空间
	int *Center_XV = new int[ImgWidth];// 椭圆中心点x坐标参数空间
	int *Center_YV = new int[ImgHeight];// 椭圆中心点y坐标参数空间
	int *A_axis_V = new int[max(ImgWidth,ImgHeight)/2];// 椭圆长轴参数空间
	int *B_axis_V = new int[max(ImgWidth,ImgHeight)/2];// 椭圆短轴参数空间
	
	memset(Center_XV,0,sizeof(int)*ImgWidth);
	memset(Center_YV,0,sizeof(int)*ImgHeight);
	memset(A_axis_V,0,sizeof(int)*max(ImgWidth,ImgHeight)/2);
	memset(B_axis_V,0,sizeof(int)*max(ImgWidth,ImgHeight)/2);

	double Theta,X_c,Y_c,A_axis,B_axis;

int loop = 1;// 成功求出参数的迭代次数
	int looptop = loop * 1;// 总的迭代次数(也就是控制计算时间的上限,以免陷入无限循环)
while(loop > 0 && looptop > 0)
{
looptop --;
	int idx;
	for (count = totalPt; count > 0; count--)// 打乱样本点排列的顺序
	{
		POINT ptrtmp;
		idx = rand() % count;
		
		ptrtmp = seq[idx];
		seq[idx] = seq[count-1];
		seq[count-1] = ptrtmp;		
	}

	double PioMatrix[5*5];
	for (i = 0; i < 5; i++)
	{
		PioMatrix[i*5] = seq[i].x * seq[i].x;
		PioMatrix[i*5 + 1] = 2 * seq[i].x * seq[i].y;
		PioMatrix[i*5 + 2] = seq[i].y * seq[i].y;
		PioMatrix[i*5 + 3] = 2 * seq[i].x;
		PioMatrix[i*5 + 4] = 2 * seq[i].y;
	}

	if (GaussJordanInv(PioMatrix,5) == false)// Gauss-Jordan求逆
		continue;
	double sum;
	for (i = 0; i < 5; i++)
	{
		sum = 0;
		for (j = 0; j < 5; j++)
		{
			sum +=  -(PioMatrix[i*5 + j]);
		}
		Para[i] = sum;
	}

	if (pow(Para[1],2) - Para[0] * Para[2] > 0)
			continue;

		if (fabs(Para[0] - Para[2]) < 1e-20)
			Theta = 1.570796326;
		else if (Para[0] > Para[2])
			Theta = 0.5 * (atan(2.0 * Para[1] / (Para[0] - Para[2])) + PI);
		else
			Theta = 0.5 * (atan(2.0 * Para[1] / (Para[0] - Para[2])));

		X_c = (4.0 * Para[1] * Para[4] - 4.0 * Para[2] * Para[3]) / (4.0 * Para[0] * Para[2] - 4.0 * Para[1] * Para[1]);
		Y_c = (4.0 * Para[1] * Para[3] - 4.0 * Para[0] * Para[4]) / (4.0 * Para[0] * Para[2] - 4.0 * Para[1] * Para[1]);
		A_axis = 2 * (Para[0] * pow(X_c,2) + Para[2] * pow(Y_c,2) + 2 * Para[1] * X_c * Y_c - 1) 
						/ (Para[0] + Para[2] - sqrt(pow(Para[0] - Para[2],2) + pow(2.0 * Para[1],2)));
		B_axis = 2 * (Para[0] * pow(X_c,2) + Para[2] * pow(Y_c,2) + 2 * Para[1] * X_c * Y_c - 1) 
						/ (Para[0] + Para[2] + sqrt(pow(Para[0] - Para[2],2) + pow(2.0 * Para[1],2)));
		
		A_axis = sqrt(A_axis);	//长轴
		B_axis = sqrt(B_axis);	//短轴

		int AngleTmp = (int)(Theta * 180 / PI + 360 + 0.5) % 360;
		Angle_V[AngleTmp]++;
		
		if (X_c < 0 || Y_c < 0 || A_axis < 0 || B_axis < 0)
			continue;
		if (X_c >= ImgWidth || Y_c >= ImgHeight || A_axis > max(ImgWidth,ImgHeight)/2 || B_axis > max(ImgWidth,ImgHeight)/2)
			continue;

		if (X_c >= 0 && X_c < ImgWidth)
			Center_XV[(int)X_c]++;
		if (Y_c >= 0 && Y_c < ImgHeight)
			Center_YV[(int)Y_c]++;
		if (A_axis >= 0 && A_axis < max(ImgWidth,ImgHeight)/2)
			A_axis_V[(int)A_axis]++;
		if (B_axis >= 0 && B_axis < max(ImgWidth,ImgHeight)/2)
			B_axis_V[(int)B_axis]++;		
		loop--;
}
	
	int Angle,Ai,Bi,Cx,Cy;
	//	Angle
	int MaxPara = 0;
	for (i = 0; i < 360; i++)
	{
		if (MaxPara < Angle_V[i])
		{
			MaxPara = Angle_V[i];
			Angle = i;
		}
	}
	//	Cy
	MaxPara = 0;
	for (i = 0; i < ImgHeight; i++)
	{
		if (MaxPara < Center_YV[i])
		{
			MaxPara = Center_YV[i];
			Cy = i;
		}
	}
	//	Cx
	MaxPara = 0;
	for (i = 0; i < ImgWidth; i++)
	{
		if (MaxPara < Center_XV[i])
		{
			MaxPara = Center_XV[i];
			Cx = i;
		}
	}
	//	Ai
	MaxPara = 0;
	for (i = 0; i < max(ImgWidth,ImgHeight)/2; i++)
	{
		if (MaxPara < A_axis_V[i])
		{
			MaxPara = A_axis_V[i];
			Ai = i;
		}
	}
	//	Bi
	MaxPara = 0;
	for (i = 0; i < max(ImgWidth,ImgHeight)/2; i++)
	{
		if (MaxPara < B_axis_V[i])
		{
			MaxPara = B_axis_V[i];
			Bi = i;
		}
	}

	delete[] Center_XV;
	delete[] Center_YV;
	delete[] A_axis_V;
	delete[] B_axis_V;


	double sma = SinMem[Angle];
	double cma = CosMem[Angle];
	for (int n = 0; n < 360; n++)
	{
		i = (Bi) * CosMem[360 - n];
		j = (Ai) * SinMem[360 - n];
		
		int x,y;
		x = (j * cma - i * sma) + Cx;
		y = (i * cma + j * sma) + Cy;
		
		Mask[y * ImgWidth + x] = 0;
	}
	delete[] Mask;
	delete[] seq;
}

测试结果:

原图:

拟合结果(虚线为拟合的椭圆):

前面说到这种方法有缺陷,请看下面的情形:

原图:

拟合结果:

       当样本点只集中在椭圆的一边时,随机5点的hough变换总会拟合错误,实际应用中往往会发生这样的情况。这是因为公式错了吗?于是我单独提取出五点做测试,即只做一次迭代。测试结果如下图所示。图中实线为实际椭圆,我是用画图工具拖出来的“完美椭圆”,用橡皮擦擦掉一大半部分,最后再做一点旋转。打交叉的是取样的5点,虚线是用这5点代入公式求得的拟合椭圆。可见求得的椭圆穿过了5个点,表明不是求解错了,可就是跟实际的有很大差别。唯一的解释是取样点不在我们想要的椭圆上,也就是说即使是用画图工具拖出来看似完美的椭圆并不完美,这是因为样本点的坐标是整型,精度很低。所以随机5点的hough变换存在很严重的系统误差,当取样点分散在椭圆上下左右时,这种误差会比较小,当集中在某个区域时,误差就会非常大。

   

        解决办法就是多采几个点,然后用最小二乘法求解。

下图是10点随机采样的结果:

下图是将所有点一起计算的结果:

可见,采样点越多,拟合度越好。但是一次取样的点越多,这些点落入相同椭圆的概率就越小。这就需要一些手段把椭圆的边缘从噪声中提取出来。至于最小二乘法的椭圆检测算法我将另开一贴来讨论。