Skip to content
Closed
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
139 changes: 98 additions & 41 deletions src/phg/matching/descriptor_matcher.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -8,14 +8,26 @@ void phg::DescriptorMatcher::filterMatchesRatioTest(const std::vector<std::vecto
{
filtered_matches.clear();

throw std::runtime_error("not implemented yet");
}
constexpr float distanceRatioThreshold = 0.75f;

for (const std::vector<cv::DMatch>& match : matches) {
if (match.size() > 1) {
if (match[0].distance / match[1].distance <= distanceRatioThreshold) {
filtered_matches.push_back(match.front());
}
} else if (match.size() == 1) {
filtered_matches.push_back(match.front());
} else {
std::cerr << "Descriptor matcher ratio test: got empty match?" << std::endl;
}
}
}

void phg::DescriptorMatcher::filterMatchesClusters(const std::vector<cv::DMatch> &matches,
const std::vector<cv::KeyPoint> keypoints_query,
const std::vector<cv::KeyPoint> keypoints_train,
std::vector<cv::DMatch> &filtered_matches)
std::vector<cv::DMatch> &filtered_matches,
int iterations)
{
filtered_matches.clear();

Expand All @@ -35,42 +47,87 @@ 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-дерево
// NOTE: Точное дерево, может, и можно построить, но похоже, что OpenCV
// в этом случае этого всё-таки не делает. Для большей точности лишь
// используем увеличенное количество проверок (на результаты тестов не особо
// влияет, но при использовании нескольких итераций кластерного фильтра
// увеличенное число проверок действительно даёт более стабильный результат).
std::shared_ptr<cv::flann::IndexParams> index_params = flannKdTreeIndexParams(4);
std::shared_ptr<cv::flann::SearchParams> search_params = flannKsTreeSearchParams(64);

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 общих элементов
for (int i_query = 0; i_query < n_matches; ++i_query) {
size_t curr_consistent_matches = 1;

for (int j_query = 1; j_query < total_neighbours; ++j_query) {
if (distances2_query.at<float>(i_query, j_query) > radius2_query) {
break;
}

// поскольку points_query и points_train заполнены по порядку соответствующими точками
const int target_point_index = indices_query.at<int>(i_query, j_query);
const int i_train = i_query;

for (int j_train = 1; j_train < total_neighbours; ++j_train) {
if (distances2_train.at<float>(i_train, j_train) > radius2_train) {
break;
}

if (indices_train.at<int>(i_train, j_train) == target_point_index) {
++curr_consistent_matches;
break;
}
}
}

if (curr_consistent_matches >= consistent_matches) {
filtered_matches.push_back(matches[i_query]);
}
}

// NOTE: реализация применения фильтра в несколько итераций — для
// минимизации количества изменений в коде задания сделано рекурсивно,
// в чистовом варианте можно (и лучше) заменить циклом.
if (iterations > 1) {
if (filtered_matches.size() >= total_neighbours) {
std::vector<cv::DMatch> maches_pass2 = filtered_matches;
filterMatchesClusters(maches_pass2, keypoints_query, keypoints_train, filtered_matches, iterations - 1);
} else {
// Слишком мало точек для следующей итерации.
filtered_matches.clear();
}
}
}
3 changes: 2 additions & 1 deletion src/phg/matching/descriptor_matcher.h
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,8 @@ namespace phg {
static void filterMatchesClusters(const std::vector<cv::DMatch> &matches,
const std::vector<cv::KeyPoint> keypoints_query,
const std::vector<cv::KeyPoint> keypoints_train,
std::vector<cv::DMatch> &filtered_matches);
std::vector<cv::DMatch> &filtered_matches,
int iterations = 3);
};

}
24 changes: 21 additions & 3 deletions src/phg/matching/flann_matcher.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,11 @@
phg::FlannMatcher::FlannMatcher()
{
// параметры для приближенного поиска
// index_params = flannKdTreeIndexParams(TODO);
// search_params = flannKsTreeSearchParams(TODO);
// NOTE: Поскольку более обоснованных идей у меня нет, здесь используются
// стандартные значения по умолчанию из OpenCV. Похоже, что они работают
// вполне неплохо.
index_params = flannKdTreeIndexParams(4);
search_params = flannKsTreeSearchParams(32);
}

void phg::FlannMatcher::train(const cv::Mat &train_desc)
Expand All @@ -17,5 +20,20 @@ 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;
cv::Mat dists;
flann_index->knnSearch(query_desc, indices, dists, k);

matches.clear();
matches.reserve(indices.rows);

for (int r = 0; r < indices.rows; ++r) {
std::vector<cv::DMatch>& match = matches.emplace_back();
match.reserve(indices.cols);

for (int c = 0; c < indices.cols; ++c) {
const float dist = std::sqrt(dists.at<float>(r, c)); // knnSearch возвращает квадрат расстояния
match.push_back(cv::DMatch(r, indices.at<int>(r, c), dist));
}
}
}
120 changes: 66 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({0, 0, 0, -x0 * w1, -y0 * w1, -w0 * w1, x0 * y1, y0 * y1, -w0 * y1});
A.push_back({x0 * w1, y0 * w1, w0 * w1, 0, 0, 0, -x0 * x1, -y0 * x1, w0 * x1});
}

int res = gauss(A, H);
Expand Down Expand Up @@ -168,57 +168,66 @@ 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
// NOTE: n_samples следует оставить равным 4, поскольку именно столько
// точек нужно (и используется ниже в коде) для оценки гомографии.
// Количество попыток 800 примерно соответствует вероятности успеха в
// 99.9% и доле корректных точек 0.3. В конечном итоге самая низкая
// доля "голосов" в тестах за лучшую гомографию была около 0.35 (727/2070),
// но при более низком количестве попыток был зафиксирован и результат,
// лишь немного превышающий 0.3 (638/2070). Поэтому с запасом можно
// оставить 800 попыток, если это не будет вредить производительности
// (а на тестах проблем с этим пока не заметно).

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🔥

const int n_trials = 800;
const int n_samples = 4;

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 +247,10 @@ 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");
const double x = T.at<double>(0, 0) * pt.x + T.at<double>(0, 1) * pt.y + T.at<double>(0, 2);
const double y = T.at<double>(1, 0) * pt.x + T.at<double>(1, 1) * pt.y + T.at<double>(1, 2);
const double w = T.at<double>(2, 0) * pt.x + T.at<double>(2, 1) * pt.y + T.at<double>(2, 2);
return cv::Point2d(x / w, y / w);
}

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

for (size_t i = 0; i < n_images; ++i) {
if (parent[i] == -1) {
cv::Mat identity_mat(3, 3, CV_64FC1, 0.0);
identity_mat.at<double>(0, 0) = 1.0;
identity_mat.at<double>(1, 1) = 1.0;
identity_mat.at<double>(2, 2) = 1.0;
parent_homographies[i] = std::move(identity_mat);
} else {
parent_homographies[i] = homography_builder(imgs[i], imgs[parent[i]]);
}
}

// NOTE: Тут можно модифицировать алгоритм, чтобы не перемножать лишний
// раз матрицы гомографий для одних и тех же пар, но пока оставил здесь
// простой вариант.
for (size_t i = 0; i < n_images; ++i) {
Hs[i] = parent_homographies[i].clone();
int parent_idx = parent[i];

while (parent_idx != -1) {
Hs[i] = parent_homographies[parent_idx] * Hs[i];
parent_idx = parent[parent_idx];
}
}
}

bbox2<double, cv::Point2d> bbox;
Expand Down
Loading
Loading