光斑计算类

#pragma once
/*--------------------------------------------------------------------------------------------------
	*  @Copyright (c) , All rights reserved.
	*  @file:       svm_train.h
	*  @version:    ver 1.0
	*  @author:     闹闹
	*  @brief:
	*  @change:
	*  @email: 	1319144981@qq.com
	*  Date             Version    Changed By      Changes
	*  2021/10/27 9:21:33    1.0       闹闹            create
	写一句自己最喜欢的话吧。
	为天地立心，为生民立命，为往圣继绝学，为万世开太平。
--------------------------------------------------------------------------------------------------*/
#include <iostream>
#include <Eigen/Dense>
#include "opencv2/opencv.hpp"
#include "opencv2/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <algorithm>
#ifndef __GAUSS_SPOT__H__
#define __GAUSS_SPOT__H__


namespace nao {
	namespace algorithm{
		namespace gauss_spot {
		
			/**
			* @brief 初始化数据矩阵
			* @param img      单通道float型数据
			* @param vector_A 向量A
			* @param matrix_B 矩阵B
			* @return
			* @参考链接 https://blog.csdn.net/houjixin/article/details/8490653
			*/
			bool initData(cv::Mat img, Eigen::VectorXf& vector_A, Eigen::MatrixXf& matrix_B) {
				if (img.empty())
					return false;
				int length = img.cols * img.rows;
				Eigen::VectorXf tmp_A(length);
				Eigen::MatrixXf tmp_B(length, 5);
				int i = 0, j = 0, iPos = 0;

				while (i < img.cols) {
					j = 0;
					while (j < img.rows) {
						float value = *img.ptr<float>(i, j);
						if (value <= 0.5) value = 1;
						tmp_A(iPos) = value * log(value);
						tmp_B(iPos, 0) = value;
						tmp_B(iPos, 1) = value * i;
						tmp_B(iPos, 2) = value * j;
						tmp_B(iPos, 3) = value * i * i;
						tmp_B(iPos, 4) = value * j * j;
						++iPos;
						++j;
					}
					++i;
				}
				vector_A = tmp_A;
				matrix_B = tmp_B;
			}
			/**
			 * @brief 获取光斑中心
			 * @param x0
			 * @param y0
			 * @param vector_A
			 * @param matrix_B
			 * @return
			*/
			bool getCentrePoint(float& x0, float& y0, const Eigen::VectorXf vector_A, const Eigen::MatrixXf matrix_B) {
				//QR分解
				Eigen::HouseholderQR<Eigen::MatrixXf> qr;
				qr.compute(matrix_B);
				Eigen::MatrixXf R = qr.matrixQR().triangularView<Eigen::Upper>();
				Eigen::MatrixXf Q = qr.householderQ();

				//块操作，获取向量或矩阵的局部
				Eigen::VectorXf S;
				S = (Q.transpose() * vector_A).head(5);
				Eigen::MatrixXf R1;
				R1 = R.block(0, 0, 5, 5);

				Eigen::VectorXf C;
				C = R1.inverse() * S;

				x0 = -0.5 * C[1] / C[3];
				y0 = -0.5 * C[2] / C[4];
				return true;
			}

			/**
			 * @brief 高斯拟合求光斑,求光斑的中心，测试函数
			 * @param src 
			*/
			void gauss_spot_center(const cv::Mat src, std::vector<cv::Point2f> & dst_pt) {
				//cv::Mat src = cv::imread("1650.bmp", 0);
				cv::Mat img = src.clone();
				cv::Mat dis_img;
				cv::cvtColor(img, dis_img, cv::COLOR_GRAY2BGR);

				cv::Mat thresh_img;

				//二值化
				double dCurdwTh = 200;
				cv::threshold(img, thresh_img, dCurdwTh, 255, cv::THRESH_BINARY);
				//cv::Mat element = cv::getStructuringElement(cv::MORPH_RECT, cv::Size(3, 3));;//X腐蚀
				//cv::Mat elementP = cv::getStructuringElement(cv::MORPH_RECT, cv::Size(5, 5));//膨胀
				//cv::erode(thresh_img, thresh_img, element);//腐蚀操作
				//cv::dilate(thresh_img, thresh_img, elementP);//膨胀操作

				std::vector<std::vector<cv::Point>> contours;
				std::vector<cv::Vec4i> hierarchy;

				cv::findContours(thresh_img, contours, hierarchy, cv::RETR_CCOMP, cv::CHAIN_APPROX_NONE);
				std::vector<cv::Point2f> pt;
				for (int i = 0; i < contours.size(); i++) {
					cv::Rect rect = cv::boundingRect(contours.at(i));
					if (rect.width < 4 || rect.width > 65 || rect.height < 4)
						continue;

					cv::Mat tmp = img(cv::Rect(rect.x - 5, rect.y - 5, rect.width + 10, rect.height + 10)).clone();
					tmp.convertTo(tmp, CV_32F);
					Eigen::VectorXf vector_A;
					Eigen::MatrixXf matrix_B;

					initData(tmp, vector_A, matrix_B);
					float x0 = 0;
					float y0 = 0;
					getCentrePoint(x0, y0, vector_A, matrix_B);
					cv::Point2f p;

					if (x0 != 0 && y0 != 0) {
						p.x = x0 + rect.tl().x - 5;
						p.y = y0 + rect.tl().y - 5;
						pt.emplace_back(p);
						cv::circle(dis_img, p, 1, cv::Scalar(0, 0, 255), -1);
					}
				}
				for (int i = 0; i < pt.size();i++) {
					dst_pt.emplace_back(pt[i]);
				}
				cv::waitKey(0);
			}
		}//namespace gauss_spot
	}//algorithm
}//namespace nao
#endif // !__GAUSS_SPOT__H__


/*----------------------------------------------------------------------------- (C) COPYRIGHT LEI *****END OF FILE------------------------------------------------------------------------------*/

#pragma once
/*--------------------------------------------------------------------------------------------------
	*  @Copyright (c) , All rights reserved.
	*  @file:       svm_train.h
	*  @version:    ver 1.0
	*  @author:     闹闹
	*  @brief:
	*  @change:
	*  @email: 	1319144981@qq.com
	*  Date             Version    Changed By      Changes
	*  2021/10/27 9:21:33    1.0       闹闹            create
	写一句自己最喜欢的话吧。
	为天地立心，为生民立命，为往圣继绝学，为万世开太平。
--------------------------------------------------------------------------------------------------*/

#include <iostream>
#include "opencv2/opencv.hpp"
#include "opencv2/highgui.hpp"

namespace nao {
	namespace algorithm {
		namespace hessen_light {

			/**
			 * @brief https://blog.csdn.net/dangkie/article/details/78996761
			*/
			void StegerLine(const cv::Mat src, std::vector<cv::Point2f>& dst_pt) {
				//cv::Mat src = cv::imread("1650.bmp", 0);
				cv::Mat img = src.clone();

				//高斯滤波
				img.convertTo(img, CV_32FC1);
				cv::GaussianBlur(img, img, cv::Size(0, 0), 2.5, 2.5);

				cv::Mat dx, dy;
				cv::Mat dxx, dyy, dxy;

				//一阶偏导数
				cv::Mat mDx, mDy;
				//二阶偏导数
				cv::Mat mDxx, mDyy, mDxy;

				mDx = (cv::Mat_<float>(1, 2) << 1, -1);//x偏导
				mDy = (cv::Mat_<float>(2, 1) << 1, -1);//y偏导
				mDxx = (cv::Mat_<float>(1, 3) << 1, -2, 1);//二阶x偏导
				mDyy = (cv::Mat_<float>(3, 1) << 1, -2, 1);//二阶y偏导
				mDxy = (cv::Mat_<float>(2, 2) << 1, -1, -1, 1);//二阶xy偏导

				cv::filter2D(img, dx, CV_32FC1, mDx);
				cv::filter2D(img, dy, CV_32FC1, mDy);
				cv::filter2D(img, dxx, CV_32FC1, mDxx);
				cv::filter2D(img, dyy, CV_32FC1, mDyy);
				cv::filter2D(img, dxy, CV_32FC1, mDxy);

				//hessian矩阵
				int cols = src.cols;
				int rows = src.rows;
				std::vector<cv::Point2f> pts;

				for (int col = 0; col < cols; ++col)
				{
					for (int row = rows - 1; row != -1; --row)
					{
						if (src.at<uchar>(row, col) < 210) continue;

						cv::Mat hessian(2, 2, CV_32FC1);
						hessian.at<float>(0, 0) = dxx.at<float>(row, col);
						hessian.at<float>(0, 1) = dxy.at<float>(row, col);
						hessian.at<float>(1, 0) = dxy.at<float>(row, col);
						hessian.at<float>(1, 1) = dyy.at<float>(row, col);
						cv::Mat eValue;
						cv::Mat eVectors;
						cv::eigen(hessian, eValue, eVectors);
						double nx, ny;
						double fmaxD = 0;

						if (std::fabs(eValue.at<float>(0, 0)) >= std::fabs(eValue.at<float>(1, 0)))  //求特征值最大时对应的特征向量
						{
							nx = eVectors.at<float>(0, 0);
							ny = eVectors.at<float>(0, 1);
							fmaxD = eValue.at<float>(0, 0);
						}
						else
						{
							nx = eVectors.at<float>(1, 0);
							ny = eVectors.at<float>(1, 1);
							fmaxD = eValue.at<float>(1, 0);
						}

						float t = -(nx * dx.at<float>(row, col) + ny * dy.at<float>(row, col))
							/ (nx * nx * dxx.at<float>(row, col) + 2 * nx * ny * dxy.at<float>(row, col) + ny * ny * dyy.at<float>(row, col));
						float tnx = t * nx;
						float tny = t * ny;

						if (std::fabs(tnx) <= 0.25 && std::fabs(tny) <= 0.25)
						{
							float x = col + /*.5*/tnx;
							float y = row + /*.5*/tny;
							pts.push_back({ x, y });
							break;
						}
					}
				}

				cv::Mat display;
				cv::cvtColor(src, display, cv::COLOR_GRAY2BGR);

				for (int k = 0; k < pts.size(); k++)
				{
					cv::Point rpt;
					rpt.x = pts[k].x;
					rpt.y = pts[k].y;
					cv::circle(display, rpt, 1, cv::Scalar(0, 0, 255));
					dst_pt.emplace_back(rpt);
					
				}
				cv::imshow("result", display);
				cv::waitKey(0);
			}

		}//namespace hessen_light
	}//namespace algorithm
}//namespace nao

/*----------------------------------------------------------------------------- (C) COPYRIGHT LEI *****END OF FILE------------------------------------------------------------------------------*/