|
|
@@ -0,0 +1,405 @@
|
|
|
+#include "vis_marker_detector.h"
|
|
|
+#include "cuda_helper.hpp"
|
|
|
+#include "image_process/vis_marker_kernels.cuh"
|
|
|
+#include "simple_mq.h"
|
|
|
+#include "utility.hpp"
|
|
|
+#include "variable_defs.h"
|
|
|
+
|
|
|
+#include <opencv2/calib3d.hpp>
|
|
|
+#include <opencv2/imgproc.hpp>
|
|
|
+
|
|
|
+#include <thread>
|
|
|
+
|
|
|
+using namespace simple_mq_singleton;
|
|
|
+
|
|
|
+namespace vis_marker_detector_impl {
|
|
|
+
|
|
|
+ double cross_2d(const Eigen::Vector2d &a, const Eigen::Vector2d &b) {
|
|
|
+ return a.x() * b.y() - b.x() * a.y();
|
|
|
+ }
|
|
|
+
|
|
|
+ template<typename T>
|
|
|
+ int8_t sgn(T val) {
|
|
|
+ return (T(0) < val) - (val < T(0));
|
|
|
+ }
|
|
|
+
|
|
|
+ Eigen::Vector3f to_homo_vec(const Eigen::Vector2f &vec) {
|
|
|
+ return {vec.x(), vec.y(), 1.f};
|
|
|
+ }
|
|
|
+
|
|
|
+ Eigen::Matrix3f to_intrinsic_matrix(const camera_intrinsic &cam) {
|
|
|
+ Eigen::Matrix3f ret = Eigen::Matrix3f::Zero();
|
|
|
+ ret(0, 0) = cam.fx;
|
|
|
+ ret(0, 2) = cam.cx;
|
|
|
+ ret(1, 1) = cam.fy;
|
|
|
+ ret(1, 2) = cam.cy;
|
|
|
+ ret(2, 2) = 1.f;
|
|
|
+ return ret;
|
|
|
+ }
|
|
|
+
|
|
|
+ Eigen::Matrix3f to_cross_matrix(const Eigen::Vector3f &vec) {
|
|
|
+ Eigen::Matrix3f ret;
|
|
|
+ ret << 0, -vec.z(), vec.y(),
|
|
|
+ vec.z(), 0, -vec.x(),
|
|
|
+ -vec.y(), vec.x(), 0;
|
|
|
+ return ret;
|
|
|
+ }
|
|
|
+
|
|
|
+ Eigen::Matrix3f to_essential_matrix(const Eigen::Isometry3f &transform) {
|
|
|
+ Eigen::Vector3f trans_vec = transform.translation();
|
|
|
+ Eigen::Matrix3f rot_mat = transform.rotation();
|
|
|
+ // E = np.cross([T], R)
|
|
|
+ return to_cross_matrix(trans_vec) * rot_mat;
|
|
|
+ }
|
|
|
+
|
|
|
+ Eigen::Matrix3f to_fundamental_matrix(const stereo_camera_info &info) {
|
|
|
+ auto essential_mat = to_essential_matrix(info.transform);
|
|
|
+ auto left_mat = to_intrinsic_matrix(*info.left);
|
|
|
+ auto right_mat = to_intrinsic_matrix(*info.right);
|
|
|
+ // F = np.linalg.inv(np.transpose(K2)) @ E @ np.linalg.inv(K1)
|
|
|
+ return right_mat.transpose().inverse()
|
|
|
+ * essential_mat
|
|
|
+ * left_mat.inverse();
|
|
|
+ }
|
|
|
+
|
|
|
+ void undistort_points(const camera_intrinsic &cam, Eigen::Matrix2Xf *points) {
|
|
|
+ Eigen::Matrix3f in_mat = to_intrinsic_matrix(cam);
|
|
|
+ Eigen::Matrix3f in_mat_inv = in_mat.inverse();
|
|
|
+ for (int i = 0; i < points->cols(); ++i) {
|
|
|
+ auto p = points->col(i);
|
|
|
+ auto p_n = Eigen::Vector2f{ // p_n = np.linalg.inv(K) @ p
|
|
|
+ (p.x() - cam.cx) / cam.fx,
|
|
|
+ (p.y() - cam.cy) / cam.fy};
|
|
|
+ auto r2 = p_n.x() * p_n.x() + p_n.y() * p_n.y();
|
|
|
+ auto factor = 1 + cam.k[0] * r2 + cam.k[1] * r2 * r2;
|
|
|
+ p_n *= factor;
|
|
|
+ p = Eigen::Vector2f{
|
|
|
+ p_n.x() * cam.fx + cam.cx,
|
|
|
+ p_n.y() * cam.fy + cam.cy};
|
|
|
+ (*points)(0, i) = p.x();
|
|
|
+ (*points)(1, i) = p.y();
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ // compatible with OpenCV
|
|
|
+ using project_matrix_type = Eigen::Matrix<float, 3, 4, Eigen::RowMajor>;
|
|
|
+
|
|
|
+ project_matrix_type to_proj_mat(const camera_intrinsic &cam,
|
|
|
+ const Eigen::Isometry3f &transform) {
|
|
|
+ return to_intrinsic_matrix(cam) * transform.matrix().block<3, 4>(0, 0);
|
|
|
+ }
|
|
|
+
|
|
|
+ // undistorted points in pixel
|
|
|
+ void match_points(const Eigen::Matrix2Xf &points_left,
|
|
|
+ const Eigen::Matrix2Xf &points_right,
|
|
|
+ const stereo_camera_info &cam,
|
|
|
+ Eigen::Matrix3Xf *out,
|
|
|
+ float max_tolerance = 1.f) {
|
|
|
+ static constexpr auto max_points = 10;
|
|
|
+ static constexpr auto reserve_size = 2 * max_points;
|
|
|
+
|
|
|
+ assert(points_left.cols() <= max_points);
|
|
|
+ assert(points_right.cols() <= max_points);
|
|
|
+
|
|
|
+ auto match_left = std::vector<cv::Point2f>{};
|
|
|
+ auto match_right = std::vector<cv::Point2f>{};
|
|
|
+ match_left.reserve(reserve_size);
|
|
|
+ match_right.reserve(reserve_size);
|
|
|
+
|
|
|
+ auto f_mat = to_fundamental_matrix(cam);
|
|
|
+ for (auto lp: points_left.colwise())
|
|
|
+ for (auto rp: points_right.colwise()) {
|
|
|
+ auto p_l = to_homo_vec(lp);
|
|
|
+ auto p_r = to_homo_vec(rp);
|
|
|
+ auto dis = p_l.transpose() * f_mat * p_r;
|
|
|
+ if (abs(dis.value()) < max_tolerance) {
|
|
|
+ match_left.emplace_back(lp.x(), lp.y());
|
|
|
+ match_right.emplace_back(rp.x(), rp.y());
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ auto match_cnt = match_left.size();
|
|
|
+ auto result = Eigen::MatrixX4f{match_cnt, 4};
|
|
|
+ auto ret_mat = cv::Mat{4, (int) match_cnt, CV_32FC1, result.data()};
|
|
|
+ auto proj_left = to_proj_mat(*cam.left, Eigen::Isometry3f::Identity());
|
|
|
+ auto proj_right = to_proj_mat(*cam.right, cam.transform);
|
|
|
+ auto proj_left_mat = cv::Mat{3, 4, CV_32FC1, proj_left.data()};
|
|
|
+ auto proj_right_mat = cv::Mat{3, 4, CV_32FC1, proj_right.data()};
|
|
|
+ cv::triangulatePoints(proj_left_mat, proj_right_mat,
|
|
|
+ match_left, match_right,
|
|
|
+ ret_mat);
|
|
|
+
|
|
|
+ out->resize(3, match_cnt);
|
|
|
+ for (int k = 0; k < match_cnt; ++k)
|
|
|
+ for (int i = 0; i < 3; ++i) { // homo to non-homo
|
|
|
+ (*out)(i, k) = result(k, i) / result(k, 3);
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ float calc_transform_error(const Eigen::Matrix3Xf &points,
|
|
|
+ const Eigen::Matrix3Xf &marker,
|
|
|
+ const Eigen::Isometry3f &trans) {
|
|
|
+ Eigen::VectorXf err{marker.cols()};
|
|
|
+ for (auto k = 0; k < marker.cols(); ++k) {
|
|
|
+ auto mp = marker.col(k);
|
|
|
+ auto tp = trans * mp;
|
|
|
+ err(k) = std::numeric_limits<float>::max();
|
|
|
+ for (auto p: points.colwise()) {
|
|
|
+ auto dis = (p - tp).norm();
|
|
|
+ err(k) = std::min(dis, err(k));
|
|
|
+ }
|
|
|
+ }
|
|
|
+// return err.maxCoeff();
|
|
|
+ return err.norm() / sqrtf(err.size()); // RMS
|
|
|
+ }
|
|
|
+
|
|
|
+ std::optional<Eigen::Isometry3f>
|
|
|
+ match_marker(const Eigen::Matrix3Xf &points, const Eigen::Matrix3Xf &marker,
|
|
|
+ float *err = nullptr, float max_tolerance = 1.f) {
|
|
|
+ assert(marker.cols() >= 3);
|
|
|
+ if (points.cols() < marker.cols()) return {};
|
|
|
+
|
|
|
+ float dis1 = (marker.col(1) - marker.col(0)).norm();
|
|
|
+ float dis2 = (marker.col(2) - marker.col(0)).norm();
|
|
|
+ float dis3 = (marker.col(2) - marker.col(1)).norm();
|
|
|
+ auto m_rep = marker.block<3, 3>(0, 0);
|
|
|
+ for (auto p0: points.colwise())
|
|
|
+ for (auto p1: points.colwise()) {
|
|
|
+ if (abs((p1 - p0).norm() - dis1) > max_tolerance) [[likely]] continue;
|
|
|
+ for (auto p2: points.colwise()) {
|
|
|
+ if (abs((p2 - p0).norm() - dis2) > max_tolerance) [[likely]] continue;
|
|
|
+ if (abs((p2 - p1).norm() - dis3) > max_tolerance) [[likely]] continue;
|
|
|
+ Eigen::Matrix3f p_rep;
|
|
|
+ p_rep << p0, p1, p2;
|
|
|
+ auto trans = Eigen::Isometry3f{Eigen::umeyama(m_rep, p_rep)};
|
|
|
+ if (err != nullptr) {
|
|
|
+ *err = calc_transform_error(points, marker, trans);
|
|
|
+ }
|
|
|
+ return trans;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return {};
|
|
|
+ }
|
|
|
+
|
|
|
+ struct smart_detect_kernels {
|
|
|
+ smart_gpu_buffer<int8_t> kernel_dev;
|
|
|
+ smart_gpu_buffer<uint16_t> kernel_weight_dev;
|
|
|
+
|
|
|
+ void create(uint8_t border, uint8_t kernel_num) {
|
|
|
+ if (border == last_border && kernel_num == last_num) return;
|
|
|
+ generate_kernels(border, kernel_num);
|
|
|
+ }
|
|
|
+
|
|
|
+ kernel_bunch<int8_t> to_kernel_bunch() {
|
|
|
+ return {last_border, last_num, kernel_dev.ptr};
|
|
|
+ }
|
|
|
+
|
|
|
+ private:
|
|
|
+ smart_buffer<uint16_t> kernel_weight;
|
|
|
+ uint8_t last_border = 0, last_num = 0;
|
|
|
+
|
|
|
+ void generate_kernels(uint8_t border, uint8_t kernel_num) {
|
|
|
+ int width = 2 * border + 1;
|
|
|
+ int size[] = {width, width, kernel_num};
|
|
|
+ auto ker_mat = cv::Mat{3, size, CV_8SC1};
|
|
|
+ kernel_weight.create(kernel_num);
|
|
|
+ for (int k = 0; k < kernel_num; ++k) {
|
|
|
+ auto th = std::numbers::pi / 2 / kernel_num * k;
|
|
|
+ auto vec0 = Eigen::Vector2d{std::cos(th), std::sin(th)};
|
|
|
+ auto vec1 = Eigen::Vector2d{-std::sin(th), std::cos(th)};
|
|
|
+ kernel_weight.ptr[k] = 0;
|
|
|
+ for (int i = 0; i < width; ++i) {
|
|
|
+ int x = i - border;
|
|
|
+ for (int j = 0; j < width; ++j) {
|
|
|
+ int y = j - border;
|
|
|
+ auto vec = Eigen::Vector2d{x, y};
|
|
|
+ auto cross1 = cross_2d(vec, vec0);
|
|
|
+ auto cross2 = cross_2d(vec, vec1);
|
|
|
+ if (std::abs(cross1) < 1 || std::abs(cross2) < 1) {
|
|
|
+ ker_mat.at<int8_t>(j, i, k) = 0;
|
|
|
+ } else {
|
|
|
+ auto val = sgn(cross1 * cross2);
|
|
|
+ ker_mat.at<int8_t>(j, i, k) = val;
|
|
|
+ if (val > 0) {
|
|
|
+ ++kernel_weight.ptr[k];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ kernel_dev.upload_from(ker_mat.data, ker_mat.total());
|
|
|
+ kernel_weight_dev.upload_from(kernel_weight);
|
|
|
+ last_border = border;
|
|
|
+ last_num = kernel_num;
|
|
|
+ }
|
|
|
+ };
|
|
|
+
|
|
|
+ struct corner_detector {
|
|
|
+
|
|
|
+ smart_detect_kernels *kernels;
|
|
|
+ smart_cuda_stream stream;
|
|
|
+
|
|
|
+ cv::cuda::GpuMat score_dev;
|
|
|
+ smart_gpu_buffer<ushort2> corner_coord;
|
|
|
+ smart_gpu_buffer<uint32_t> corner_index;
|
|
|
+ smart_buffer<ushort2> corner_coord_host;
|
|
|
+
|
|
|
+ struct options {
|
|
|
+ float score_threshold = 0.2;
|
|
|
+ uint16_t max_corners = 10;
|
|
|
+ uint8_t cuda_block_width = 16;
|
|
|
+ };
|
|
|
+
|
|
|
+ // accept 8-bit (CV_8UC1) or 12-bit (CV_16UC1)
|
|
|
+ void detect(const cv::cuda::GpuMat &in, const cv::Mat &in_host,
|
|
|
+ Eigen::Matrix2Xf *out, const options &opt) {
|
|
|
+ constexpr auto block_width = 16;
|
|
|
+ constexpr auto lm_border = 1;
|
|
|
+ auto refine_criteria = cv::TermCriteria{
|
|
|
+ cv::TermCriteria::COUNT | cv::TermCriteria::EPS, 10, 0.001};
|
|
|
+ auto refine_win_size = cv::Size{5, 5};
|
|
|
+ auto refine_zero_zone = cv::Size{-1, -1};
|
|
|
+
|
|
|
+ assert(in.channels() == 1);
|
|
|
+ assert(in_host.channels() == 1);
|
|
|
+
|
|
|
+ score_dev.create(in.size(), CV_8UC1);
|
|
|
+ auto score_img = to_image_type<uint8_t>(score_dev);
|
|
|
+
|
|
|
+ // score evaluate
|
|
|
+ switch (in.depth()) {
|
|
|
+ case CV_8U: {
|
|
|
+ auto in_img = to_image_type<uint8_t>(in);
|
|
|
+ call_corner_score(in_img, score_img,
|
|
|
+ kernels->to_kernel_bunch(),
|
|
|
+ kernels->kernel_weight_dev.ptr,
|
|
|
+ block_width, stream.cuda);
|
|
|
+ break;
|
|
|
+ }
|
|
|
+ case CV_16U: {
|
|
|
+ auto in_img = to_image_type<uint16_t>(in);
|
|
|
+ call_corner_score(in_img, score_img,
|
|
|
+ kernels->to_kernel_bunch(),
|
|
|
+ kernels->kernel_weight_dev.ptr,
|
|
|
+ block_width, stream.cuda);
|
|
|
+ break;
|
|
|
+ }
|
|
|
+ default: {
|
|
|
+ RET_ERROR;
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ // corner filter
|
|
|
+ corner_coord.create(opt.max_corners);
|
|
|
+ corner_index.create(1);
|
|
|
+ static constexpr uint32_t zero = 0;
|
|
|
+ corner_index.upload_from(&zero, sizeof(uint32_t), stream.cuda); // corner_index = 0;
|
|
|
+ auto filter_result = filter_result_type<ushort2>{
|
|
|
+ opt.max_corners, corner_index.ptr, corner_coord.ptr};
|
|
|
+ call_lm_filter(score_img, filter_result,
|
|
|
+ opt.score_threshold, lm_border, block_width, stream.cuda);
|
|
|
+ uint32_t corner_count = 0;
|
|
|
+ CUDA_API_CHECK(cudaMemcpyAsync(&corner_count, corner_index.ptr,
|
|
|
+ sizeof(uint32_t), cudaMemcpyDeviceToHost, stream.cuda));
|
|
|
+ corner_coord.download_to(&corner_coord_host, stream.cuda);
|
|
|
+ CUDA_API_CHECK(cudaStreamSynchronize(stream.cuda));
|
|
|
+ if (corner_count > opt.max_corners) {
|
|
|
+ corner_count = opt.max_corners;
|
|
|
+ SPDLOG_WARN("Try increase #corner limit.");
|
|
|
+ }
|
|
|
+ out->resize(2, corner_count);
|
|
|
+ for (int i = 0; i < corner_count; ++i) {
|
|
|
+ (*out)(0, i) = corner_coord_host.ptr[i].x;
|
|
|
+ (*out)(1, i) = corner_coord_host.ptr[i].y;
|
|
|
+ }
|
|
|
+ if (corner_count == 0) return; // no corner find
|
|
|
+
|
|
|
+ // corner refinement
|
|
|
+ auto coord_mat = cv::Mat{(int) corner_count, 2, CV_32FC1, out->data()};
|
|
|
+ cv::Mat in_mat;
|
|
|
+ switch (in_host.depth()) {
|
|
|
+ case CV_8U: {
|
|
|
+ in_mat = in_host;
|
|
|
+ break;
|
|
|
+ }
|
|
|
+ case CV_16U: {
|
|
|
+ constexpr auto factor = 1.f / (1 << 4); // assume 12-bit to 8-bit conversion
|
|
|
+ in_host.convertTo(in_mat, CV_8U, factor);
|
|
|
+ break;
|
|
|
+ }
|
|
|
+ default: {
|
|
|
+ RET_ERROR;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ cv::cornerSubPix(in_mat, coord_mat,
|
|
|
+ refine_win_size, refine_zero_zone, refine_criteria);
|
|
|
+ }
|
|
|
+
|
|
|
+ };
|
|
|
+
|
|
|
+}
|
|
|
+
|
|
|
+using namespace vis_marker_detector_impl;
|
|
|
+
|
|
|
+struct vis_marker_detector::impl {
|
|
|
+
|
|
|
+ smart_detect_kernels kernels;
|
|
|
+ corner_detector detector_left, detector_right;
|
|
|
+ Eigen::Matrix2Xf corner_left, corner_right;
|
|
|
+ Eigen::Matrix3Xf space_points;
|
|
|
+
|
|
|
+ std::optional<Eigen::Isometry3f>
|
|
|
+ detect(smart_mat *left_img, smart_mat *right_img,
|
|
|
+ const stereo_camera_info &cam_info, const Eigen::Matrix3Xf &marker,
|
|
|
+ const options &opt, float *err = nullptr) {
|
|
|
+
|
|
|
+ kernels.create(opt.detect_kernel_border,
|
|
|
+ opt.detect_kernel_number);
|
|
|
+ detector_left.kernels = &kernels;
|
|
|
+ detector_right.kernels = &kernels;
|
|
|
+
|
|
|
+ corner_detector::options corner_opt;
|
|
|
+ corner_opt.max_corners = 2 * marker.cols();
|
|
|
+ corner_opt.score_threshold = opt.score_threshold;
|
|
|
+ corner_opt.cuda_block_width = opt.cuda_block_width;
|
|
|
+
|
|
|
+ // detect corners in multiple threads
|
|
|
+ auto help_t = std::thread([&] {
|
|
|
+ // binding cuda context
|
|
|
+ auto cuda_ctx = mq().query_variable<CUcontext>(CUDA_CONTEXT);
|
|
|
+ CUDA_API_CHECK(cuCtxSetCurrent(cuda_ctx));
|
|
|
+ detector_left.detect(left_img->dev(detector_left.stream.cv),
|
|
|
+ left_img->host(detector_left.stream.cv),
|
|
|
+ &corner_left, corner_opt);
|
|
|
+ });
|
|
|
+ detector_right.detect(right_img->dev(detector_right.stream.cv),
|
|
|
+ right_img->host(detector_right.stream.cv),
|
|
|
+ &corner_right, corner_opt);
|
|
|
+ help_t.join();
|
|
|
+
|
|
|
+ // no sufficient corners for this marker
|
|
|
+ if (corner_left.cols() < marker.cols()
|
|
|
+ || corner_right.cols() < marker.cols()) {
|
|
|
+ return {};
|
|
|
+ }
|
|
|
+
|
|
|
+ undistort_points(*cam_info.left, &corner_left);
|
|
|
+ undistort_points(*cam_info.right, &corner_right);
|
|
|
+ match_points(corner_left, corner_right, cam_info, &space_points, opt.pixel_tolerance);
|
|
|
+ return match_marker(space_points, marker, err, opt.space_tolerance);
|
|
|
+ }
|
|
|
+
|
|
|
+};
|
|
|
+
|
|
|
+vis_marker_detector::vis_marker_detector()
|
|
|
+ : pimpl(std::make_unique<impl>()) {}
|
|
|
+
|
|
|
+vis_marker_detector::~vis_marker_detector() = default;
|
|
|
+
|
|
|
+std::optional<Eigen::Isometry3f>
|
|
|
+vis_marker_detector::detect(smart_mat *left_img, smart_mat *right_img,
|
|
|
+ const stereo_camera_info &cam_info, const Eigen::Matrix3Xf &marker,
|
|
|
+ const vis_marker_detector::options &opt, float *err) {
|
|
|
+ return pimpl->detect(left_img, right_img, cam_info, marker, opt, err);
|
|
|
+}
|