Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
119 changes: 76 additions & 43 deletions src/phg/matching/descriptor_matcher.cpp
Original file line number Diff line number Diff line change
@@ -1,17 +1,21 @@
#include "descriptor_matcher.h"

#include <opencv2/flann/miniflann.hpp>
#include "flann_factory.h"
#include <opencv2/flann/miniflann.hpp>
#include <unordered_set>

void phg::DescriptorMatcher::filterMatchesRatioTest(const std::vector<std::vector<cv::DMatch>> &matches,
std::vector<cv::DMatch> &filtered_matches)
{
filtered_matches.clear();

throw std::runtime_error("not implemented yet");
const float threshold = 0.6;
for (auto match : matches) {
if (match[0].distance < threshold * match[1].distance) {
filtered_matches.push_back(match[0]);
}
}
}


void phg::DescriptorMatcher::filterMatchesClusters(const std::vector<cv::DMatch> &matches,
const std::vector<cv::KeyPoint> keypoints_query,
const std::vector<cv::KeyPoint> keypoints_train,
Expand All @@ -35,42 +39,71 @@ void phg::DescriptorMatcher::filterMatchesClusters(const std::vector<cv::DMatch>
points_query.at<cv::Point2f>(i) = keypoints_query[matches[i].queryIdx].pt;
points_train.at<cv::Point2f>(i) = keypoints_train[matches[i].trainIdx].pt;
}
//
// // размерность всего 2, так что точное KD-дерево
// std::shared_ptr<cv::flann::IndexParams> index_params = flannKdTreeIndexParams(TODO);
// std::shared_ptr<cv::flann::SearchParams> search_params = flannKsTreeSearchParams(TODO);
//
// std::shared_ptr<cv::flann::Index> index_query = flannKdTreeIndex(points_query, index_params);
// std::shared_ptr<cv::flann::Index> index_train = flannKdTreeIndex(points_train, index_params);
//
// // для каждой точки найти total neighbors ближайших соседей
// cv::Mat indices_query(n_matches, total_neighbours, CV_32SC1);
// cv::Mat distances2_query(n_matches, total_neighbours, CV_32FC1);
// cv::Mat indices_train(n_matches, total_neighbours, CV_32SC1);
// cv::Mat distances2_train(n_matches, total_neighbours, CV_32FC1);
//
// index_query->knnSearch(points_query, indices_query, distances2_query, total_neighbours, *search_params);
// index_train->knnSearch(points_train, indices_train, distances2_train, total_neighbours, *search_params);
//
// // оценить радиус поиска для каждой картинки
// // NB: radius2_query, radius2_train: квадраты радиуса!
// float radius2_query, radius2_train;
// {
// std::vector<double> max_dists2_query(n_matches);
// std::vector<double> max_dists2_train(n_matches);
// for (int i = 0; i < n_matches; ++i) {
// max_dists2_query[i] = distances2_query.at<float>(i, total_neighbours - 1);
// max_dists2_train[i] = distances2_train.at<float>(i, total_neighbours - 1);
// }
//
// int median_pos = n_matches / 2;
// std::nth_element(max_dists2_query.begin(), max_dists2_query.begin() + median_pos, max_dists2_query.end());
// std::nth_element(max_dists2_train.begin(), max_dists2_train.begin() + median_pos, max_dists2_train.end());
//
// radius2_query = max_dists2_query[median_pos] * radius_limit_scale * radius_limit_scale;
// radius2_train = max_dists2_train[median_pos] * radius_limit_scale * radius_limit_scale;
// }
//
// метч остается, если левое и правое множества первых total_neighbors соседей в радиусах поиска(radius2_query, radius2_train) имеют как минимум consistent_matches общих элементов
// // TODO заполнить filtered_matches
}

// размерность всего 2, так что точное KD-дерево
std::shared_ptr<cv::flann::IndexParams> index_params = flannKdTreeIndexParams(4);
std::shared_ptr<cv::flann::SearchParams> search_params = flannKsTreeSearchParams(32);

std::shared_ptr<cv::flann::Index> index_query = flannKdTreeIndex(points_query, index_params);
std::shared_ptr<cv::flann::Index> index_train = flannKdTreeIndex(points_train, index_params);

// для каждой точки найти total neighbors ближайших соседей
cv::Mat indices_query(n_matches, total_neighbours, CV_32SC1);
cv::Mat distances2_query(n_matches, total_neighbours, CV_32FC1);
cv::Mat indices_train(n_matches, total_neighbours, CV_32SC1);
cv::Mat distances2_train(n_matches, total_neighbours, CV_32FC1);

index_query->knnSearch(points_query, indices_query, distances2_query, total_neighbours, *search_params);
index_train->knnSearch(points_train, indices_train, distances2_train, total_neighbours, *search_params);

// оценить радиус поиска для каждой картинки
// NB: radius2_query, radius2_train: квадраты радиуса!
float radius2_query, radius2_train;
{
std::vector<double> max_dists2_query(n_matches);
std::vector<double> max_dists2_train(n_matches);
for (int i = 0; i < n_matches; ++i) {
max_dists2_query[i] = distances2_query.at<float>(i, total_neighbours - 1);
max_dists2_train[i] = distances2_train.at<float>(i, total_neighbours - 1);
}

int median_pos = n_matches / 2;
std::nth_element(max_dists2_query.begin(), max_dists2_query.begin() + median_pos, max_dists2_query.end());
std::nth_element(max_dists2_train.begin(), max_dists2_train.begin() + median_pos, max_dists2_train.end());

radius2_query = max_dists2_query[median_pos] * radius_limit_scale * radius_limit_scale;
radius2_train = max_dists2_train[median_pos] * radius_limit_scale * radius_limit_scale;
}

// метч остается, если левое и правое множества первых total_neighbours соседей в радиусах поиска(radius2_query, radius2_train)
// имеют как минимум consistent_matches общих элементов

for (int i = 0; i < n_matches; i++) {
std::unordered_set<int> in_radius_query;
int match_count = 0;

// находим все индексы в радиусе radius2_query
for (int j = 0; j < total_neighbours; j++) {
int ind_query = indices_query.at<int>(i, j);
float dist2_query = distances2_query.at<float>(i, j);
if (dist2_query <= radius2_query) {
in_radius_query.insert(ind_query);
}
}

// считаем индексы в радиусе radius2_train, существующие в in_radius_query
for (int j = 0; j < total_neighbours; j++) {
int ind_train = indices_query.at<int>(i, j);
float dist2_train = distances2_train.at<float>(i, j);
bool exists_in_query = in_radius_query.find(ind_train) != in_radius_query.end();
if (dist2_train <= radius2_train && exists_in_query) {
match_count++;
}
}

// если матчей достаточно - добавляем
if (match_count >= consistent_matches) {
filtered_matches.push_back(matches[i]);
}
}
}
16 changes: 13 additions & 3 deletions src/phg/matching/flann_matcher.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,8 @@
phg::FlannMatcher::FlannMatcher()
{
// параметры для приближенного поиска
// index_params = flannKdTreeIndexParams(TODO);
// search_params = flannKsTreeSearchParams(TODO);
index_params = flannKdTreeIndexParams(4);
search_params = flannKsTreeSearchParams(32);
}

void phg::FlannMatcher::train(const cv::Mat &train_desc)
Expand All @@ -17,5 +17,15 @@ void phg::FlannMatcher::train(const cv::Mat &train_desc)

void phg::FlannMatcher::knnMatch(const cv::Mat &query_desc, std::vector<std::vector<cv::DMatch>> &matches, int k) const
{
throw std::runtime_error("not implemented yet");
cv::Mat indices, dists;
flann_index->knnSearch(query_desc, indices, dists, k, *search_params);
matches.resize(query_desc.rows);

for (int i = 0; i < query_desc.rows; i++) {
for (int j = 0; j < k; j++) {
int index = indices.at<int>(i, j);
float dist = dists.at<float>(i, j);
matches[i].emplace_back(i, index, dist);
}
}
}
118 changes: 64 additions & 54 deletions src/phg/sfm/homography.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -84,8 +84,8 @@ namespace {
double w1 = ws1[i];

// 8 elements of matrix + free term as needed by gauss routine
// A.push_back({TODO});
// A.push_back({TODO});
A.push_back({ x0, y0, w0, 0, 0, 0, -x1 * x0, -x1 * y0, x1 * w0 });
A.push_back({ 0, 0, 0, x0, y0, w0, -y1 * x0, -y1 * y0, y1 * w0 });
}

int res = gauss(A, H);
Expand Down Expand Up @@ -168,57 +168,57 @@ namespace {
// * (простое описание для понимания)
// * [3] http://ikrisoft.blogspot.com/2015/01/ransac-with-contrario-approach.html

// const int n_matches = points_lhs.size();
//
// // https://en.wikipedia.org/wiki/Random_sample_consensus#Parameters
// const int n_trials = TODO;
//
// const int n_samples = TODO;
// uint64_t seed = 1;
// const double reprojection_error_threshold_px = 2;
//
// int best_support = 0;
// cv::Mat best_H;
//
// std::vector<int> sample;
// for (int i_trial = 0; i_trial < n_trials; ++i_trial) {
// randomSample(sample, n_matches, n_samples, &seed);
//
// cv::Mat H = estimateHomography4Points(points_lhs[sample[0]], points_lhs[sample[1]], points_lhs[sample[2]], points_lhs[sample[3]],
// points_rhs[sample[0]], points_rhs[sample[1]], points_rhs[sample[2]], points_rhs[sample[3]]);
//
// int support = 0;
// for (int i_point = 0; i_point < n_matches; ++i_point) {
// try {
// cv::Point2d proj = phg::transformPoint(points_lhs[i_point], H);
// if (cv::norm(proj - cv::Point2d(points_rhs[i_point])) < reprojection_error_threshold_px) {
// ++support;
// }
// } catch (const std::exception &e)
// {
// std::cerr << e.what() << std::endl;
// }
// }
//
// if (support > best_support) {
// best_support = support;
// best_H = H;
//
// std::cout << "estimateHomographyRANSAC : support: " << best_support << "/" << n_matches << std::endl;
//
// if (best_support == n_matches) {
// break;
// }
// }
// }
//
// std::cout << "estimateHomographyRANSAC : best support: " << best_support << "/" << n_matches << std::endl;
//
// if (best_support == 0) {
// throw std::runtime_error("estimateHomographyRANSAC : failed to estimate homography");
// }
//
// return best_H;
const int n_matches = points_lhs.size();

// https://en.wikipedia.org/wiki/Random_sample_consensus#Parameters
const int n_trials = 1000;

const int n_samples = 10;
uint64_t seed = 1;
const double reprojection_error_threshold_px = 2;

int best_support = 0;
cv::Mat best_H;

std::vector<int> sample;
for (int i_trial = 0; i_trial < n_trials; ++i_trial) {
randomSample(sample, n_matches, n_samples, &seed);

cv::Mat H = estimateHomography4Points(points_lhs[sample[0]], points_lhs[sample[1]], points_lhs[sample[2]], points_lhs[sample[3]],
points_rhs[sample[0]], points_rhs[sample[1]], points_rhs[sample[2]], points_rhs[sample[3]]);

int support = 0;
for (int i_point = 0; i_point < n_matches; ++i_point) {
try {
cv::Point2d proj = phg::transformPoint(points_lhs[i_point], H);
if (cv::norm(proj - cv::Point2d(points_rhs[i_point])) < reprojection_error_threshold_px) {
++support;
}
} catch (const std::exception &e)
{
std::cerr << e.what() << std::endl;
}
}

if (support > best_support) {
best_support = support;
best_H = H;

std::cout << "estimateHomographyRANSAC : support: " << best_support << "/" << n_matches << std::endl;

if (best_support == n_matches) {
break;
}
}
}

std::cout << "estimateHomographyRANSAC : best support: " << best_support << "/" << n_matches << std::endl;

if (best_support == 0) {
throw std::runtime_error("estimateHomographyRANSAC : failed to estimate homography");
}

return best_H;
}

}
Expand All @@ -238,7 +238,17 @@ cv::Mat phg::findHomographyCV(const std::vector<cv::Point2f> &points_lhs, const
// таким преобразованием внутри занимается функции cv::perspectiveTransform и cv::warpPerspective
cv::Point2d phg::transformPoint(const cv::Point2d &pt, const cv::Mat &T)
{
throw std::runtime_error("not implemented yet");
// [x'] [T00 T01 T02] [x]
// [y'] = [T10 T11 T12] [y]
// [w'] [T20 T21 T22] [1]
double x = T.at<double>(0, 0) * pt.x + T.at<double>(0, 1) * pt.y + T.at<double>(0, 2);
double y = T.at<double>(1, 0) * pt.x + T.at<double>(1, 1) * pt.y + T.at<double>(1, 2);
double w = T.at<double>(2, 0) * pt.x + T.at<double>(2, 1) * pt.y + T.at<double>(2, 2);

double x_transformed = x / w;
double y_transformed = y / w;

return { x_transformed, y_transformed };
}

cv::Point2d phg::transformPointCV(const cv::Point2d &pt, const cv::Mat &T) {
Expand Down
56 changes: 42 additions & 14 deletions src/phg/sfm/panorama_stitcher.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,35 @@ cv::Mat phg::stitchPanorama(const std::vector<cv::Mat> &imgs,
{
// здесь надо посчитать вектор Hs
// при этом можно обойтись n_images - 1 вызовами функтора homography_builder
throw std::runtime_error("not implemented yet");

int root = -1;
for (int i = 0; i < n_images; i++) {
if (parent[i] == -1) {
root = i;
break;
}
}

Hs.resize(n_images);
std::vector<bool> done(n_images, false);

Hs[root] = cv::Mat::eye(3, 3, CV_64FC1);
done[root] = true;

std::function<void(size_t)> compute_homography = [&](size_t idx) {
if (done[idx])
return;
compute_homography(parent[idx]);
cv::Mat cur_img = imgs[idx];
cv::Mat parent_img = imgs[parent[idx]];
cv::Mat homography = homography_builder(cur_img, parent_img);
Hs[idx] = Hs[parent[idx]] * homography;
done[idx] = true;
};

for (int i = 0; i < n_images; i++) {
compute_homography(i);
}
}

bbox2<double, cv::Point2d> bbox;
Expand All @@ -46,19 +74,19 @@ cv::Mat phg::stitchPanorama(const std::vector<cv::Mat> &imgs,
// из-за растяжения пикселей при использовании прямой матрицы гомографии после отображения между пикселями остается пустое пространство
// лучше использовать обратную и для каждого пикселя на итоговвой картинке проверять, с какой картинки он может получить цвет
// тогда в некоторых пикселях цвет будет дублироваться, но изображение будет непрерывным
// for (int i = 0; i < n_images; ++i) {
// for (int y = 0; y < imgs[i].rows; ++y) {
// for (int x = 0; x < imgs[i].cols; ++x) {
// cv::Vec3b color = imgs[i].at<cv::Vec3b>(y, x);
//
// cv::Point2d pt_dst = applyH(cv::Point2d(x, y), Hs[i]) - bbox.min();
// int y_dst = std::max(0, std::min((int) std::round(pt_dst.y), result_height - 1));
// int x_dst = std::max(0, std::min((int) std::round(pt_dst.x), result_width - 1));
//
// result.at<cv::Vec3b>(y_dst, x_dst) = color;
// }
// }
// }
// for (int i = 0; i < n_images; ++i) {
// for (int y = 0; y < imgs[i].rows; ++y) {
// for (int x = 0; x < imgs[i].cols; ++x) {
// cv::Vec3b color = imgs[i].at<cv::Vec3b>(y, x);

// cv::Point2d pt_dst = applyH(cv::Point2d(x, y), Hs[i]) - bbox.min();
// int y_dst = std::max(0, std::min((int) std::round(pt_dst.y), result_height - 1));
// int x_dst = std::max(0, std::min((int) std::round(pt_dst.x), result_width - 1));

// result.at<cv::Vec3b>(y_dst, x_dst) = color;
// }
// }
// }

std::vector<cv::Mat> Hs_inv;
std::transform(Hs.begin(), Hs.end(), std::back_inserter(Hs_inv), [&](const cv::Mat &H){ return H.inv(); });
Expand Down
Loading
Loading