From a93773a8a352eb02f00f1a755da7c3326964afdc Mon Sep 17 00:00:00 2001 From: Ekaterina Date: Fri, 20 Mar 2026 23:58:26 +0300 Subject: [PATCH 1/2] task02 --- src/phg/matching/descriptor_matcher.cpp | 119 +++++++++++++++--------- src/phg/matching/flann_matcher.cpp | 16 +++- src/phg/sfm/homography.cpp | 118 ++++++++++++----------- src/phg/sfm/panorama_stitcher.cpp | 56 ++++++++--- 4 files changed, 195 insertions(+), 114 deletions(-) diff --git a/src/phg/matching/descriptor_matcher.cpp b/src/phg/matching/descriptor_matcher.cpp index f4bcd87..b6c4ae3 100644 --- a/src/phg/matching/descriptor_matcher.cpp +++ b/src/phg/matching/descriptor_matcher.cpp @@ -1,17 +1,21 @@ #include "descriptor_matcher.h" -#include #include "flann_factory.h" +#include +#include void phg::DescriptorMatcher::filterMatchesRatioTest(const std::vector> &matches, std::vector &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 &matches, const std::vector keypoints_query, const std::vector keypoints_train, @@ -35,42 +39,71 @@ void phg::DescriptorMatcher::filterMatchesClusters(const std::vector points_query.at(i) = keypoints_query[matches[i].queryIdx].pt; points_train.at(i) = keypoints_train[matches[i].trainIdx].pt; } -// -// // размерность всего 2, так что точное KD-дерево -// std::shared_ptr index_params = flannKdTreeIndexParams(TODO); -// std::shared_ptr search_params = flannKsTreeSearchParams(TODO); -// -// std::shared_ptr index_query = flannKdTreeIndex(points_query, index_params); -// std::shared_ptr 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 max_dists2_query(n_matches); -// std::vector max_dists2_train(n_matches); -// for (int i = 0; i < n_matches; ++i) { -// max_dists2_query[i] = distances2_query.at(i, total_neighbours - 1); -// max_dists2_train[i] = distances2_train.at(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 index_params = flannKdTreeIndexParams(4); + std::shared_ptr search_params = flannKsTreeSearchParams(32); + + std::shared_ptr index_query = flannKdTreeIndex(points_query, index_params); + std::shared_ptr 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 max_dists2_query(n_matches); + std::vector max_dists2_train(n_matches); + for (int i = 0; i < n_matches; ++i) { + max_dists2_query[i] = distances2_query.at(i, total_neighbours - 1); + max_dists2_train[i] = distances2_train.at(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 in_radius_query; + int match_count = 0; + + // находим все индексы в радиусе radius2_query + for (int j = 0; j < total_neighbours; j++) { + int ind_query = indices_query.at(i, j); + float dist2_query = distances2_query.at(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(i, j); + float dist2_train = distances2_train.at(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]); + } + } +} \ No newline at end of file diff --git a/src/phg/matching/flann_matcher.cpp b/src/phg/matching/flann_matcher.cpp index 9e9f518..d692b5a 100644 --- a/src/phg/matching/flann_matcher.cpp +++ b/src/phg/matching/flann_matcher.cpp @@ -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) @@ -17,5 +17,15 @@ void phg::FlannMatcher::train(const cv::Mat &train_desc) void phg::FlannMatcher::knnMatch(const cv::Mat &query_desc, std::vector> &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(i, j); + float dist = dists.at(i, j); + matches[i].emplace_back(i, index, dist); + } + } } diff --git a/src/phg/sfm/homography.cpp b/src/phg/sfm/homography.cpp index 5cbc780..aa8fcd0 100644 --- a/src/phg/sfm/homography.cpp +++ b/src/phg/sfm/homography.cpp @@ -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); @@ -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 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 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; } } @@ -238,7 +238,17 @@ cv::Mat phg::findHomographyCV(const std::vector &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(0, 0) * pt.x + T.at(0, 1) * pt.y + T.at(0, 2); + double y = T.at(1, 0) * pt.x + T.at(1, 1) * pt.y + T.at(1, 2); + double w = T.at(2, 0) * pt.x + T.at(2, 1) * pt.y + T.at(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) { diff --git a/src/phg/sfm/panorama_stitcher.cpp b/src/phg/sfm/panorama_stitcher.cpp index 8d76939..586bfb9 100644 --- a/src/phg/sfm/panorama_stitcher.cpp +++ b/src/phg/sfm/panorama_stitcher.cpp @@ -23,7 +23,35 @@ cv::Mat phg::stitchPanorama(const std::vector &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); + bool done[n_images] = {}; + + Hs[root] = cv::Mat::eye(3, 3, CV_64FC1); + done[root] = true; + + std::function 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 bbox; @@ -46,19 +74,19 @@ cv::Mat phg::stitchPanorama(const std::vector &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(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(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(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(y_dst, x_dst) = color; + // } + // } + // } std::vector Hs_inv; std::transform(Hs.begin(), Hs.end(), std::back_inserter(Hs_inv), [&](const cv::Mat &H){ return H.inv(); }); From 000984543290092d0a9f42f6a9d9a5e0e5c479c1 Mon Sep 17 00:00:00 2001 From: Ekaterina Date: Sat, 21 Mar 2026 00:07:21 +0300 Subject: [PATCH 2/2] task02 --- src/phg/sfm/panorama_stitcher.cpp | 2 +- src/phg/sift/sift.cpp | 295 ++++++++++++++++++------------ tests/test_matching.cpp | 2 +- tests/test_sift.cpp | 26 +-- 4 files changed, 197 insertions(+), 128 deletions(-) diff --git a/src/phg/sfm/panorama_stitcher.cpp b/src/phg/sfm/panorama_stitcher.cpp index 586bfb9..93e8727 100644 --- a/src/phg/sfm/panorama_stitcher.cpp +++ b/src/phg/sfm/panorama_stitcher.cpp @@ -33,7 +33,7 @@ cv::Mat phg::stitchPanorama(const std::vector &imgs, } Hs.resize(n_images); - bool done[n_images] = {}; + std::vector done(n_images, false); Hs[root] = cv::Mat::eye(3, 3, CV_64FC1); done[root] = true; diff --git a/src/phg/sift/sift.cpp b/src/phg/sift/sift.cpp index 7204771..9db8ce0 100755 --- a/src/phg/sift/sift.cpp +++ b/src/phg/sift/sift.cpp @@ -76,26 +76,44 @@ cv::Mat phg::toGray32F(const cv::Mat& img) return gray_float; } +// как видно по количеству комментариев, осознавались октавы долго и печально :D std::vector phg::buildOctaves(const cv::Mat& img, const phg::SIFTParams& p, int verbose_level) { - const int s = p.n_octave_layers; - const double sigma0 = p.sigma; - // взятое с потолка значение блюра который уже есть в картинке. используем для того, чтобы не так сильно блюрить базовую картинку и не терять лишний раз фичи - // upd: хотя llm не соглашается со "взятое с потолка": - // It is strictly not taken from the ceiling. sigma=0.5 is the theoretical minimum blur needed to prevent aliasing (Nyquist frequency) - // when sampling a continuous signal into a discrete grid. It is a mathematically grounded assumption for digital images. + // взятое с потолка значение блюра который уже есть в картинке + // используем для того, чтобы не так сильно блюрить базовую картинку и не терять лишний раз фичи // общая идея в том, что у нас есть какой-то сигнал реального мира, и есть входное изображение // сигнал реального мира: потенциально высочайшего разрешения, можем зумиться почти до молекул, сигма почти нулевая - // сигнал с камеры, входное изображение: было произведено усреднение хотя бы по отдельным пикселям матрицы камеры, что соответствует сигме в полпикселя + // сигнал с камеры, входное изображение: было произведено усреднение хотя бы по отдельным пикселям матрицы камеры, + // что соответствует сигме в полпикселя + + // 0.5 тут как будто бы "реальное" значение + // если мы ставим значение больше, мы говорим, что текущая картинка уже сильно заблюрена, поэтому для соответствия sigma0 будем проводить условный апскейл const double sigma_nominal = p.upscale_first ? 1.0 /*2x от неапскейленного*/ : 0.5; - const int n_layers = s + 3; // нужно +2 слоя для того чтобы крайних было по соседу для поиска максимума в scale space, и еще +1 слой, чтобы получить s DoG слоев (DoG = разность двух) - int n_octaves = std::max(1, (int)std::round(std::log2(std::min(img.cols, img.rows))) - 3); // не даунскейлим дальше размера картинки в 16 пикселей, там уже не имеет смысла что-то детектировать + // количество слоев в октаве по настройке + const int s = p.n_octave_layers; + + // размытие, которое должно быть у первого изображения + const double sigma0 = p.sigma; + + // нужно +2 слоя для того чтобы крайних было по соседу для поиска максимума в scale space, и еще +1 слой, чтобы получить s DoG слоев (DoG = разность двух) + const int n_layers = s + 3; + + // количество октав в зависимости от размера изображения + // не даунскейлим дальше размера картинки в 16 пикселей, там уже не имеет смысла что-то детектировать + int n_octaves = std::max(1, (int)std::round(std::log2(std::min(img.cols, img.rows))) - 3); + // база для слоя октавы - результат прошлого слоя cv::Mat base; + + // сигма - параметр для размытия с учетом изначального размытия изображения + // sigma_total^2 = sigma1^2 + sigma2^2 double sigma_base = std::sqrt(sigma0 * sigma0 - sigma_nominal * sigma_nominal); // можно использовать дальше как идею для инкрементального блюра слоев + + // cv::Size() - "CV, вычисли размер ядра автоматически на основе сигмы" cv::GaussianBlur(img, base, cv::Size(), sigma_base, sigma_base); + // тут храним все октавы std::vector octaves(n_octaves); for (int o = 0; o < n_octaves; o++) { @@ -104,20 +122,23 @@ std::vector phg::buildOctaves(const cv::Mat& img, const phg:: oct.layers[0] = base.clone(); - // для простоты в каждой октаве будем каждый раз блюрить базовую картинку с полной сигмой + // для простоты в каждой октаве будем каждый раз блюрить базовую картинку с полной сигмой // можно подумать, как сделать эффективнее - для построения n+1 слоя доблюревать уже поблюренный n-ый слой, так чтобы в итоге получилась такая же сигма // это будет немного быстрее, тк нужно более маленькое ядро свертки на каждый шаг for (int i = 1; i < n_layers; i++) { - // TODO double sigma_layer = sigma0 * корень из двух нужной степени, чтобы при i==s получали удвоение базового блюра; - // // вычтем sigma0 чтобы размыть ровно до нужной суммарной сигмы - // TODO sigma_layer = ... (вычитаем как в sigma base); - // cv::GaussianBlur(oct.layers[0], oct.layers[i], cv::Size(), sigma_layer, sigma_layer); + // double sigma_layer = sigma0 * корень из двух нужной степени, чтобы при i==s получали удвоение базового блюра; + // вычтем sigma0 чтобы размыть ровно до нужной суммарной сигмы + double k = std::pow(2.0, (double)i / s); // корень из двух нужной степени + double sigma_layer = sigma0 * k; + sigma_layer = std::sqrt(sigma_layer * sigma_layer - sigma0 * sigma0); + cv::GaussianBlur(oct.layers[0], oct.layers[i], cv::Size(), sigma_layer, sigma_layer); } // подготавливаем базовый слой для следующей октавы if (o + 1 < n_octaves) { // используется в opencv, формула для пересчета ключевых точек: pt_upscaled = 2^o * pt_downscaled - // TODO cv::resize(даунскейлим текущий слой в два раза, без интерполяции, просто сабсепмлинг); + // cv::resize(даунскейлим текущий слой в два раза, без интерполяции, просто сабсепмлинг); + cv::resize(oct.layers[s], base, cv::Size(), 0.5, 0.5, cv::INTER_NEAREST); // можно использовать и downsample2x_avg(oct.layers[s]), это позволяет потом заапскейлить слои обратно до оригинального разрешения без сдвига // но потребуется везде изменить формулу для пересчета ключевых точек: pt_upscaled = (pt_downscaled + 0.5) * 2^o - 0.5 @@ -136,9 +157,13 @@ std::vector phg::buildDoG(const std::vector phg::buildDoG(const std::vector phg::findScaleSpaceExtrema(const std::vector& dog, const phg::SIFTParams& params, int verbose_level) { + // количество слоев в октаве const int s = params.n_octave_layers; const double sigma0 = params.sigma; const double contrast_threshold = params.contrast_threshold; @@ -156,12 +182,12 @@ std::vector phg::findScaleSpaceExtrema(const std::vector keypoints; - for (int o = 0; o < (int)dog.size(); o++) { int real_octave = o + first_octave; @@ -193,7 +219,7 @@ std::vector phg::findScaleSpaceExtrema(const std::vector phg::findScaleSpaceExtrema(const std::vector phg::findScaleSpaceExtrema(const std::vector(yi, xi + 1) + cL.at(yi, xi - 1) - 2.f * resp_center; -// float dyy = TODO; -// float dss = TODO; -// -// float dxy = (cL.at(yi + 1, xi + 1) - cL.at(yi + 1, xi - 1) - cL.at(yi - 1, xi + 1) + cL.at(yi - 1, xi - 1)) * 0.25f; -// float dxs = TODO; -// float dys = TODO; + dxx = cL.at(yi, xi + 1) + cL.at(yi, xi - 1) - 2.f * resp_center; + dyy = cL.at(yi + 1, xi) + cL.at(yi - 1, xi) - 2.f * resp_center; + dss = pL.at(yi, xi) + nL.at(yi, xi) - 2.f * resp_center; + + dxy = (cL.at(yi + 1, xi + 1) - cL.at(yi + 1, xi - 1) - cL.at(yi - 1, xi + 1) + cL.at(yi - 1, xi - 1)) * 0.25f; + dxs = (nL.at(yi, xi + 1) - nL.at(yi, xi - 1) - pL.at(yi, xi + 1) + pL.at(yi, xi - 1)) * 0.25f; + dys = (nL.at(yi + 1, xi) - nL.at(yi - 1, xi) - pL.at(yi + 1, xi) + pL.at(yi - 1, xi)) * 0.25f; cv::Matx33f H(dxx, dxy, dxs, dxy, dyy, dys, dxs, dys, dss); cv::Vec3f g(dx, dy, ds); - // в нашей точке производная (градиент) еще не равна нулю (т.к. еще мы скорее всего не точно в оптимуме) + // в нашей точке производная (градиент) еще не равен нулю (т.к. еще мы скорее всего не точно в оптимуме) // хотим найти такой offset, где ноль производной. в предположении что оптимизируемая функция это парабола, ищем корни ее производной, линейной функции // grad(x + offset) = grad(x) + grad'(x) * offset = grad(x) + hessian(x) * offset = 0 // hessian(x) * offset = -grad(x) // линейная система. можно решить специализированным решателем либо просто найти обратную матрицу гессиана и домножить на минус градиент // offset = -hessian(x)^-1 * grad(x) cv::Vec3f offset; + // выходим, если решения нет if (!cv::solve(H, -g, offset, cv::DECOMP_LU)) break; + // ура! экстремум находится в пределах текущего пикселя if (std::abs(offset[0]) < 0.5f && std::abs(offset[1]) < 0.5f && std::abs(offset[2]) < 0.5f) { // функцию респонза оптимизировали как параболу: D(x + offset) = D(x) + grad(x) * offset + 1/2 * offset_transposed * hessian(x) * offset // подставляем hessian(x) * offset = -grad(x): D(x + offset) = D(x) + grad(x) * offset - 1/2 * offset_transposed * grad = D(x) + 1/2 * grad(x) * offset float response_optimized = resp_center + 0.5f * g.dot(offset); + + // выходим, если после уточнения так и не получили достаточно "выдающееся" значение в DoG if (std::abs(response_optimized) < thresh) break; @@ -272,22 +326,25 @@ std::vector phg::findScaleSpaceExtrema(const std::vector= ((r + 1.0) * (r + 1.0) / r)) + break; } // скейлим координаты точек обратно до родных размеров картинки @@ -323,6 +380,7 @@ std::vector phg::findScaleSpaceExtrema(const std::vector s || xi < border || xi >= cols - border || yi < border || yi >= rows - border) break; } @@ -347,13 +405,13 @@ std::vector phg::computeOrientations(const std::vector histogram(n_bins); - std::vector oriented_kpts; const int first_octave = params.upscale_first ? -1 : 0; for (const cv::KeyPoint& kp : kpts) { + std::vector histogram(n_bins); + int layer = kp.class_id; int real_octave = kp.octave; int o = real_octave - first_octave; // индекс в массиве octaves @@ -366,7 +424,8 @@ std::vector phg::computeOrientations(const std::vector phg::computeOrientations(const std::vector(py, px + 1) - img.at(py, px - 1); -// float gy = img.at(py + 1, px) - img.at(py - 1, px); -// -// float mag = TODO; -// float angle = std::atan2(TODO); // [-pi, pi] -// -// float angle_deg = angle * 180.f / (float) CV_PI; -// if (angle_deg < 0.f) angle_deg += 360.f; -// -// // гауссово взвешивание голоса точки с затуханием к краям -// float weight = std::exp(-(TODO) / (2.f * sigma_win * sigma_win)); -// if (!params.enable_orientation_gaussian_weighting) { -// weight = 1.f; -// } -// -// // голосуем в гистограмме направлений. находим два ближайших бина и гладко распределяем голос между ними -// // в таком случае, голос попавший близко к границе между бинами, проголосует поровну за оба бина -// float bin = TODO; -// if (bin >= n_bins) bin -= n_bins; -// int bin0 = (int) bin; -// int bin1 = (bin0 + 1) % n_bins; -// -// float frac = bin - bin0; -// if (!params.enable_orientation_bin_interpolation) { -// frac = 0.f; -// } -// -// histogram[bin0] += TODO; -// histogram[bin1] += TODO; + int px = xi + dx; + int py = yi + dy; + + // градиент + float gx = img.at(py, px + 1) - img.at(py, px - 1); + float gy = img.at(py + 1, px) - img.at(py - 1, px); + + float mag = std::sqrt(gx * gx + gy * gy); + float angle = std::atan2(gy, gx); // [-pi, pi] + + float angle_deg = angle * 180.f / (float)CV_PI; + if (angle_deg < 0.f) + angle_deg += 360.f; + + // гауссово взвешивание голоса точки с затуханием к краям + float weight = std::exp(-(dy * dy + dx * dx) / (2.f * sigma_win * sigma_win)); + if (!params.enable_orientation_gaussian_weighting) { + weight = 1.f; + } + + // голосуем в гистограмме направлений. находим два ближайших бина и гладко распределяем голос между ними + // в таком случае, голос попавший близко к границе между бинами, проголосует поровну за оба бина + float bin = angle_deg / (360.0 / n_bins); + if (bin >= n_bins) + bin -= n_bins; + int bin0 = (int)bin; + int bin1 = (bin0 + 1) % n_bins; + + // было ощущение, что что-то здесь не так, и мы даем бОльший вес менее приоритетному углу... :) + float frac = bin - bin0; + if (!params.enable_orientation_bin_interpolation) { + frac = 1.f; + } + + histogram[bin0] += weight * mag * frac; + histogram[bin1] += weight * mag * (1.0 - frac); } } @@ -428,7 +490,7 @@ std::vector phg::computeOrientations(const std::vector phg::computeOrientations(const std::vector a = (left + right - 2 * center) / 2 // f(1) - f(-1) = 2b -> b = (right - left) / 2 -// float offset = TODO; -// if (!params.enable_orientation_subpixel_localization) { -// offset = 0.f; -// } -// -// float bin_real = i + offset; -// if (bin_real < 0.f) bin_real += n_bins; -// if (bin_real >= n_bins) bin_real -= n_bins; -// -// float angle = bin_real * 360.f / n_bins; -// -// cv::KeyPoint new_kp = kp; -// new_kp.angle = angle; -// oriented_kpts.push_back(new_kp); + float numerator = (left - right) / 2.0; + float denominator = left + right - 2.0 * center; + + float offset = numerator / denominator; + if (!params.enable_orientation_subpixel_localization) { + offset = 0.f; + } + + float bin_real = i + offset; + if (bin_real < 0.f) + bin_real += n_bins; + if (bin_real >= n_bins) + bin_real -= n_bins; + + float angle = bin_real * 360.f / n_bins; + + cv::KeyPoint new_kp = kp; + new_kp.angle = angle; + oriented_kpts.push_back(new_kp); } } } @@ -474,7 +541,8 @@ std::vector phg::computeOrientations(const std::vector> phg::computeDescriptors(const std::vector& kpts, const std::vector& octaves, const phg::SIFTParams& params, int verbose_level) { const int s = params.n_octave_layers; @@ -486,7 +554,8 @@ std::pair> phg::computeDescriptors(const std: const int n_orient_bins = 8; const int n_dims = n_spatial_bins * n_spatial_bins * n_orient_bins; // 128 - // размер одной клетки патча в сигмах. всего размер контекста для одного дексриптора = n_spatial_bins * spatial_bin_width_sigmas сигм + // размер одной клетки патча в сигмах + // всего размер контекста для одного дексриптора = n_spatial_bins * spatial_bin_width_sigmas сигм const float spatial_bin_width_sigmas = 3.f; // в сигмах const float mag_cap = 0.2f; @@ -574,11 +643,11 @@ std::pair> phg::computeDescriptors(const std: bin_o -= n_orient_bins; // семплы вблизи края патча взвешиваем с меньшим весом -// float weight = std::exp(-(TODO) / (2.f * sigma_desc * sigma_desc)); -// if (!params.enable_descriptor_gaussian_weighting) { -// weight = 1.f; -// } -// float weighted_mag = mag * weight; + float weight = std::exp(-(rot_x * rot_x + rot_y * rot_y) / (2.f * sigma_desc * sigma_desc)); + if (!params.enable_descriptor_gaussian_weighting) { + weight = 1.f; + } + float weighted_mag = mag * weight; if (params.enable_descriptor_bin_interpolation) { // размажем вклад weighted_mag по пространственным бинам и по бинам гистограммок трилинейной интерполяцией @@ -609,8 +678,8 @@ std::pair> phg::computeDescriptors(const std: io += n_orient_bins; float wo = (dio == 0) ? (1.f - fo) : fo; -// int idx = TODO; -// desc[idx] += TODO; + int idx = (iy * n_spatial_bins + ix) * n_orient_bins + io; + desc[idx] += weighted_mag * wx * wy * wo; } } } @@ -620,9 +689,8 @@ std::pair> phg::computeDescriptors(const std: int io_nearest = (int)std::round(bin_o) % n_orient_bins; if (ix_nearest >= 0 && ix_nearest < n_spatial_bins && iy_nearest >= 0 && iy_nearest < n_spatial_bins) { - // TODO uncomment -// int idx = (iy_nearest * n_spatial_bins + ix_nearest) * n_orient_bins + io_nearest; -// desc[idx] += weighted_mag; + int idx = (iy_nearest * n_spatial_bins + ix_nearest) * n_orient_bins + io_nearest; + desc[idx] += weighted_mag; } } } @@ -712,6 +780,7 @@ void phg::SIFT::detectAndCompute(const cv::Mat& img, const cv::Mat& mask, std::v savePyramid("pyramidDoG/03_dog_octave", dog, true); kpts = findScaleSpaceExtrema(dog, p, verbose_level); + // ориентация ключевых точек это довольно дорогая операция // в случае если пользователь просит малое количество лучших точек (например, 1000, а без порога нашлось 20000), // то по производительности очень оправдано сразу их здесь и выбрать, чтобы не тащить до самого конца где все равно отбросим @@ -782,4 +851,4 @@ void phg::SIFT::savePyramid(const std::string& name, const std::vector& saveImg(ss.str(), layer); } } -} +} \ No newline at end of file diff --git a/tests/test_matching.cpp b/tests/test_matching.cpp index adaac65..214a89e 100644 --- a/tests/test_matching.cpp +++ b/tests/test_matching.cpp @@ -19,7 +19,7 @@ // TODO enable both toggles for testing custom detector & matcher -#define ENABLE_MY_DESCRIPTOR 0 +#define ENABLE_MY_DESCRIPTOR 1 #define ENABLE_MY_MATCHING 0 #define ENABLE_GPU_BRUTEFORCE_MATCHER 0 diff --git a/tests/test_sift.cpp b/tests/test_sift.cpp index cf3bd7d..54f3ece 100755 --- a/tests/test_sift.cpp +++ b/tests/test_sift.cpp @@ -28,7 +28,7 @@ // TODO ENABLE ME // TODO ENABLE ME // TODO ENABLE ME -#define ENABLE_MY_SIFT_TESTING 0 +#define ENABLE_MY_SIFT_TESTING 1 #define DENY_CREATE_REF_DATA 1 @@ -40,11 +40,11 @@ template MatchingPairData evaluateMatching(SIFT& sift, const cv: { std::vector kpts1; cv::Mat desc1; - sift.detectAndCompute(img1, { }, kpts1, desc1); + sift.detectAndCompute(img1, {}, kpts1, desc1); std::vector kpts2; cv::Mat desc2; - sift.detectAndCompute(img2, { }, kpts2, desc2); + sift.detectAndCompute(img2, {}, kpts2, desc2); // Brute-force matching with ratio test cv::BFMatcher matcher(cv::NORM_L2); @@ -225,10 +225,10 @@ void evaluateDetection(const cv::Mat& M, double minRecall, cv::Mat img0 = cv::Ma std::vector is_not_matched0(kps0.size(), true); // для каждой исходной точки хотим понять сопоставилась ли она std::vector is_not_matched1(kps1.size(), true); // для каждой точки с результирующей картинки хотим понять сопоставился ли с ней хоть кто-то - // эта прагма - способ распараллелить цикл на все ядра процессора (см. OpenMP parallel for) - // reduction позволяет сказать OpenMP что нужно провести редукцию суммированием для каждой из переменных: error_sum, n_matched, n_in_bounds, ... - // мы ведь хотим найти сумму по всем потокам - #pragma omp parallel for reduction(+ : error_sum, n_matched, n_in_bounds, size_ratio_sum, angle_diff_sum, desc_dist_sum, desc_rand_dist_sum) +// эта прагма - способ распараллелить цикл на все ядра процессора (см. OpenMP parallel for) +// reduction позволяет сказать OpenMP что нужно провести редукцию суммированием для каждой из переменных: error_sum, n_matched, n_in_bounds, ... +// мы ведь хотим найти сумму по всем потокам +#pragma omp parallel for reduction(+ : error_sum, n_matched, n_in_bounds, size_ratio_sum, angle_diff_sum, desc_dist_sum, desc_rand_dist_sum) for (ptrdiff_t i = 0; i < kps0.size(); ++i) { cv::Point2f p01 = ps01[i]; // взяли ожидаемую координату куда должна была перейти точка if (p01.x > 0 && p01.x < width && p01.y > 0 && p01.y < height) { @@ -247,8 +247,8 @@ void evaluateDetection(const cv::Mat& M, double minRecall, cv::Mat img0 = cv::Ma } } if (closest_j != -1 && min_error <= MAX_ACCEPTED_PIXEL_ERROR * width) { - // мы нашли что-то достаточно близкое - успех! - #pragma omp critical +// мы нашли что-то достаточно близкое - успех! +#pragma omp critical { is_not_matched0[i] = false; is_not_matched1[closest_j] = false; @@ -280,7 +280,7 @@ void evaluateDetection(const cv::Mat& M, double minRecall, cv::Mat img0 = cv::Ma // между точками в пространстве высокой размерности: #if 0 if (i % 100 == 0) { - #pragma omp critical +#pragma omp critical { std::cout << "d0: " << d0 << std::endl; std::cout << "d1: " << d1 << std::endl; @@ -568,7 +568,7 @@ std::string compareMats(const cv::Mat& a, const cv::Mat& b, const std::string& l return ss.str(); } if (a.empty() && b.empty()) - return { }; + return {}; // Convert to float64 for comparison cv::Mat af, bf; @@ -588,7 +588,7 @@ std::string compareMats(const cv::Mat& a, const cv::Mat& b, const std::string& l } } } - return { }; + return {}; } // ── Octave serialization ────────────────────────────────────────────────── @@ -937,4 +937,4 @@ TEST(SIFT, PairMatching) std::cout << "ENABLE_MY_SIFT_TESTING is disabled, test skipped" << std::endl; std::cout << "Final score: UNKNOWN" << std::endl; #endif -} +} \ No newline at end of file