From c53e5a7beeb756a8dd84ef5259d67da9914bb695 Mon Sep 17 00:00:00 2001 From: firelion9 Date: Sat, 21 Feb 2026 14:20:14 +0300 Subject: [PATCH 1/3] [01] SIFT impl --- .gitignore | 3 + src/phg/sift/sift.cpp | 211 ++++++++++++++++++++++-------------------- tests/test_sift.cpp | 5 +- 3 files changed, 116 insertions(+), 103 deletions(-) diff --git a/.gitignore b/.gitignore index 54b965f0..bfbcff14 100755 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,6 @@ .idea build cmake-build* +/.cache/ +/.venv/ +/.vscode diff --git a/src/phg/sift/sift.cpp b/src/phg/sift/sift.cpp index 72047711..0a7981c0 100755 --- a/src/phg/sift/sift.cpp +++ b/src/phg/sift/sift.cpp @@ -1,6 +1,7 @@ #include "sift.h" #include "libutils/rasserts.h" +#include #include #include #include @@ -108,16 +109,16 @@ std::vector phg::buildOctaves(const cv::Mat& img, const phg:: // можно подумать, как сделать эффективнее - для построения 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 * std::pow(2.0, (double) i / s); + // вычтем sigma0 чтобы размыть ровно до нужной суммарной сигмы + 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(oct.layers[s], base, oct.layers[0].size() / 2, 0, 0, cv::INTER_NEAREST); // можно использовать и downsample2x_avg(oct.layers[s]), это позволяет потом заапскейлить слои обратно до оригинального разрешения без сдвига // но потребуется везде изменить формулу для пересчета ключевых точек: pt_upscaled = (pt_downscaled + 0.5) * 2^o - 0.5 @@ -136,9 +137,13 @@ std::vector phg::buildDoG(const std::vector phg::findScaleSpaceExtrema(const std::vector phg::findScaleSpaceExtrema(const std::vector(yi, xi) - pL.at(yi, xi)) * 0.5f; // гессиан - float dxx, dxy, dyy, dxs, dys, dss; -// float dxx = cL.at(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; + float dxx = cL.at(yi, xi + 1) + cL.at(yi, xi - 1) - 2.f * resp_center; + float dyy = cL.at(yi + 1, xi) + cL.at(yi - 1, xi) - 2.f * resp_center; + float dss = nL.at(yi, xi) + pL.at(yi, xi) - 2.f * resp_center; + + 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 = (nL.at(yi, xi + 1) - nL.at(yi, xi - 1) - pL.at(yi, xi + 1) + pL.at(yi, xi - 1)) * 0.25f; + float 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); @@ -273,21 +284,22 @@ std::vector phg::findScaleSpaceExtrema(const std::vector edge_threshold) + break; } // скейлим координаты точек обратно до родных размеров картинки @@ -379,39 +391,39 @@ 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::hypot(gx, 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(-(dx * dx + dy * dy) / (2.f * sigma_win * sigma_win)); + if (!params.enable_orientation_gaussian_weighting) { + weight = 1.f; + } + + // голосуем в гистограмме направлений. находим два ближайших бина и гладко распределяем голос между ними + // в таком случае, голос попавший близко к границе между бинами, проголосует поровну за оба бина + float bin = angle_deg / 360.0f * 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 = 0.f; + } + + histogram[bin0] += (1 - frac) * weight * mag; + histogram[bin1] += frac * weight * mag; } } @@ -450,20 +462,22 @@ 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 a = (left + right - 2 * center) / 2; + float b = (right - left) / 2; + float offset = -b / (2 * a); + 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); } } } @@ -574,11 +588,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 +623,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 * wy * wx * wo; } } } @@ -620,9 +634,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; } } } diff --git a/tests/test_sift.cpp b/tests/test_sift.cpp index cf3bd7df..c660acd5 100755 --- a/tests/test_sift.cpp +++ b/tests/test_sift.cpp @@ -25,10 +25,7 @@ #define GAUSSIAN_NOISE_STDDEV 1.0 -// 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 From ae3458ab23f4d4a709fa83c892fb402561feae63 Mon Sep 17 00:00:00 2001 From: firelion9 Date: Sun, 15 Mar 2026 12:20:56 +0300 Subject: [PATCH 2/3] [02] matching base --- .gitignore | 3 +- src/phg/matching/descriptor_matcher.cpp | 147 ++++++++++++++++-------- src/phg/matching/flann_matcher.cpp | 21 +++- src/phg/sfm/homography.cpp | 129 ++++++++++++--------- src/phg/sfm/panorama_stitcher.cpp | 30 ++++- src/phg/sift/sift.cpp | 2 +- tests/test_matching.cpp | 4 +- 7 files changed, 221 insertions(+), 115 deletions(-) diff --git a/.gitignore b/.gitignore index bfbcff14..eb4e190d 100755 --- a/.gitignore +++ b/.gitignore @@ -1,7 +1,8 @@ **/*_cl.h .idea -build +/build* cmake-build* /.cache/ /.venv/ /.vscode +/data/debug/ \ No newline at end of file diff --git a/src/phg/matching/descriptor_matcher.cpp b/src/phg/matching/descriptor_matcher.cpp index f4bcd871..489f4fa7 100644 --- a/src/phg/matching/descriptor_matcher.cpp +++ b/src/phg/matching/descriptor_matcher.cpp @@ -1,27 +1,46 @@ #include "descriptor_matcher.h" -#include #include "flann_factory.h" +#include +#include +#include +#include + + /* @TUNABLE */ +static constexpr float RATIO_THRESHOLD = 0.7; +static constexpr int SEARCH_CHECKS = 10; -void phg::DescriptorMatcher::filterMatchesRatioTest(const std::vector> &matches, - std::vector &filtered_matches) +void phg::DescriptorMatcher::filterMatchesRatioTest(const std::vector>& matches, std::vector& filtered_matches) { filtered_matches.clear(); - - throw std::runtime_error("not implemented yet"); + const cv::DMatch* heap[3]; + auto cmp = [](const cv::DMatch* a, const cv::DMatch* b) { return a->distance < b->distance; }; + for (auto& match_row : matches) { + int size = 0; + for (auto& match : match_row) { + heap[size++] = &match; + std::push_heap(heap, heap + size, cmp); + if (size == 3) + std::pop_heap(heap, heap + size--, cmp); + } + if (size == 2) { + std::pop_heap(heap, heap + size, cmp); + if (heap[0]->distance / heap[1]->distance < RATIO_THRESHOLD * RATIO_THRESHOLD) { + filtered_matches.push_back(*heap[0]); + } + } else if (size == 1) { + filtered_matches.push_back(*heap[0]); + } + } } - -void phg::DescriptorMatcher::filterMatchesClusters(const std::vector &matches, - const std::vector keypoints_query, - const std::vector keypoints_train, - std::vector &filtered_matches) +void phg::DescriptorMatcher::filterMatchesClusters(const std::vector& matches, const std::vector keypoints_query, const std::vector keypoints_train, std::vector& filtered_matches) { filtered_matches.clear(); - const size_t total_neighbours = 5; // total number of neighbours to test (including candidate) - const size_t consistent_matches = 3; // minimum number of consistent matches (including candidate) - const float radius_limit_scale = 2.f; // limit search radius by scaled median + const size_t total_neighbours = 5; // total number of neighbours to test (including candidate) + const size_t consistent_matches = 3; // minimum number of consistent matches (including candidate) + const float radius_limit_scale = 2.f; // limit search radius by scaled median const int n_matches = matches.size(); @@ -35,42 +54,68 @@ 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(1); + std::shared_ptr search_params = flannKsTreeSearchParams(SEARCH_CHECKS); + + 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 общих элементов + // заполнить filtered_matches + + for (int i = 0; i < n_matches; ++i) { + int n_consistent = 0; + for (int src = 0; src < total_neighbours; ++src) { + float dsrc = distances2_query.at(i, src); + if (dsrc > radius2_query) + continue; + + for (int dst = 0; dst < total_neighbours; ++dst) { + float ddst = distances2_train.at(i, dst); + if (ddst > radius2_train) + continue; + + if (indices_query.at(i, src) == indices_train.at(i, dst)) { + n_consistent += 1; + + if (n_consistent >= consistent_matches) + goto early_success; + } + } + } + if (n_consistent >= consistent_matches) { + early_success: + filtered_matches.push_back(matches[i]); + } + } } diff --git a/src/phg/matching/flann_matcher.cpp b/src/phg/matching/flann_matcher.cpp index 9e9f5180..02c49dbc 100644 --- a/src/phg/matching/flann_matcher.cpp +++ b/src/phg/matching/flann_matcher.cpp @@ -1,13 +1,19 @@ #include +#include +#include +#include +#include #include "flann_matcher.h" #include "flann_factory.h" +static constexpr int N_TREES = 5; +static constexpr int N_CHECKS = 50; phg::FlannMatcher::FlannMatcher() { // параметры для приближенного поиска -// index_params = flannKdTreeIndexParams(TODO); -// search_params = flannKsTreeSearchParams(TODO); + index_params = flannKdTreeIndexParams(N_TREES); + search_params = flannKsTreeSearchParams(N_CHECKS); } void phg::FlannMatcher::train(const cv::Mat &train_desc) @@ -17,5 +23,14 @@ 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 idx(query_desc.rows, k, CV_32SC1); + cv::Mat dst(query_desc.rows, k, CV_32FC1); + flann_index->knnSearch(query_desc, idx, dst, k, *search_params); + matches.resize(query_desc.rows); + for (int i = 0; i < matches.size(); i++) { + matches[i].reserve(k); + for (int j = 0; j < k; ++j) { + matches[i].push_back(cv::DMatch(i, idx.at(i, j), dst.at(i, j))); + } + } } diff --git a/src/phg/sfm/homography.cpp b/src/phg/sfm/homography.cpp index 5cbc780c..7d30c770 100644 --- a/src/phg/sfm/homography.cpp +++ b/src/phg/sfm/homography.cpp @@ -1,7 +1,10 @@ #include "homography.h" +#include +#include #include #include +#include namespace { @@ -84,8 +87,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({ w1 * x0, w1 * y0, w1 * w0, 0, 0, 0, -x1 * x0, -x1 * y0, x1 * w0 }); + A.push_back({ 0, 0, 0, -w1 * x0, -w1 * y0, -w1 * w0, y1 * x0, y1 * y0, -y1 * w0 }); } int res = gauss(A, H); @@ -136,7 +139,7 @@ namespace { return *state = x; } - void randomSample(std::vector &dst, int max_id, int sample_size, uint64_t *state) + void randomSample(std::vector &dst, int max_id, int sample_size, uint64_t *state, const std::function &check_new_element) { dst.clear(); @@ -145,7 +148,7 @@ namespace { for (int i = 0; i < sample_size; ++i) { for (int k = 0; k < max_attempts; ++k) { int v = xorshift64(state) % max_id; - if (dst.empty() || std::find(dst.begin(), dst.end(), v) == dst.end()) { + if (dst.empty() || std::find(dst.begin(), dst.end(), v) == dst.end() && check_new_element(v)) { dst.push_back(v); break; } @@ -168,57 +171,68 @@ 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; + constexpr float min_dst2 = 1e-6f; + const int n_matches = points_lhs.size(); + + // https://en.wikipedia.org/wiki/Random_sample_consensus#Parameters + const int n_trials = std::log(1e-3)/std::log(1 - std::pow(0.2, 4)); + + 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 sample; + for (int i_trial = 0; i_trial < n_trials; ++i_trial) { + randomSample(sample, n_matches, n_samples, &seed, [&sample, &points_lhs, &points_rhs](int x) { + for (int prev : sample) { + auto left_vec = (points_lhs[prev] - points_lhs[x]); + if (std::hypot(left_vec.x, left_vec.y) < min_dst2) + return false; + auto right_vec = (points_rhs[prev] - points_rhs[x]); + if (std::hypot(right_vec.x, right_vec.y) < min_dst2) + return false; + } + return true; + }); + + 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 +252,10 @@ 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"); + 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); + return { x / w, y / w }; } 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 8d76939b..925b3705 100644 --- a/src/phg/sfm/panorama_stitcher.cpp +++ b/src/phg/sfm/panorama_stitcher.cpp @@ -3,6 +3,27 @@ #include #include +#include +#include + +static void callDfsOrder(const std::vector& parent, const std::function& action, std::vector& mask, int i) +{ + if (mask[i]) + return; + if (parent[i] != -1) + callDfsOrder(parent, action, mask, parent[i]); + action(i); + mask[i] = true; +} + +static void callDfsOrder(const std::vector& parent, const std::function& action) +{ + std::vector mask(parent.size()); + + for (int i = 0; i < parent.size(); ++i) { + callDfsOrder(parent, action, mask, i); + } +} /* * imgs - список картинок @@ -23,7 +44,14 @@ cv::Mat phg::stitchPanorama(const std::vector &imgs, { // здесь надо посчитать вектор Hs // при этом можно обойтись n_images - 1 вызовами функтора homography_builder - throw std::runtime_error("not implemented yet"); + callDfsOrder(parent, [&homography_builder, &parent, &imgs, &Hs, n_images](int i) { + if (parent[i] == -1) { + Hs[i] = cv::Mat({3, 3}, { 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 }); + } else { + Hs[i] = homography_builder(imgs[i], imgs[parent[i]]); + Hs[i] *= Hs[parent[i]]; + } + }); } bbox2 bbox; diff --git a/src/phg/sift/sift.cpp b/src/phg/sift/sift.cpp index 0a7981c0..0877c850 100755 --- a/src/phg/sift/sift.cpp +++ b/src/phg/sift/sift.cpp @@ -743,7 +743,7 @@ void phg::SIFT::detectAndCompute(const cv::Mat& img, const cv::Mat& mask, std::v std::tie(desc, kpts) = computeDescriptors(kpts, octaves, p, verbose_level); - // TODO всегда ли мы получаем ровно столько точек сколько запросили в параметре nfeatures? в каких случаях это не так и в какую сторону? + // всегда ли мы получаем ровно столько точек сколько запросили в параметре nfeatures? в каких случаях это не так и в какую сторону? // как подкрутить алгоритм, чтобы всегда выдавать ровно запрошенное количество точек (когда это в принципе возможно) но не сильно просесть в производительности? } diff --git a/tests/test_matching.cpp b/tests/test_matching.cpp index 158fd2c7..13657a0b 100644 --- a/tests/test_matching.cpp +++ b/tests/test_matching.cpp @@ -19,8 +19,8 @@ // TODO enable both toggles for testing custom detector & matcher -#define ENABLE_MY_DESCRIPTOR 0 -#define ENABLE_MY_MATCHING 0 +#define ENABLE_MY_DESCRIPTOR 1 +#define ENABLE_MY_MATCHING 1 #define ENABLE_GPU_BRUTEFORCE_MATCHER 0 #if ENABLE_MY_MATCHING From 1ca031b428b0a76648550765841128add1e6bacf Mon Sep 17 00:00:00 2001 From: firelion9 Date: Sun, 15 Mar 2026 13:56:04 +0300 Subject: [PATCH 3/3] [02] fixing win/macos builds --- src/phg/matching/descriptor_matcher.cpp | 2 +- src/phg/matching/flann_matcher.cpp | 4 ++-- src/phg/sfm/homography.cpp | 6 +++--- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/phg/matching/descriptor_matcher.cpp b/src/phg/matching/descriptor_matcher.cpp index 489f4fa7..f52dd8e0 100644 --- a/src/phg/matching/descriptor_matcher.cpp +++ b/src/phg/matching/descriptor_matcher.cpp @@ -7,7 +7,7 @@ #include /* @TUNABLE */ -static constexpr float RATIO_THRESHOLD = 0.7; +static constexpr float RATIO_THRESHOLD = 0.75; static constexpr int SEARCH_CHECKS = 10; void phg::DescriptorMatcher::filterMatchesRatioTest(const std::vector>& matches, std::vector& filtered_matches) diff --git a/src/phg/matching/flann_matcher.cpp b/src/phg/matching/flann_matcher.cpp index 02c49dbc..5d6d6d7a 100644 --- a/src/phg/matching/flann_matcher.cpp +++ b/src/phg/matching/flann_matcher.cpp @@ -6,8 +6,8 @@ #include "flann_matcher.h" #include "flann_factory.h" -static constexpr int N_TREES = 5; -static constexpr int N_CHECKS = 50; +static constexpr int N_TREES = 4; +static constexpr int N_CHECKS = 40; phg::FlannMatcher::FlannMatcher() { diff --git a/src/phg/sfm/homography.cpp b/src/phg/sfm/homography.cpp index 7d30c770..44be5a04 100644 --- a/src/phg/sfm/homography.cpp +++ b/src/phg/sfm/homography.cpp @@ -7,6 +7,7 @@ #include namespace { + constexpr float MIN_DST2 = 1e-6f; // источник: https://e-maxx.ru/algo/linear_systems_gauss // очень важно при выполнении метода гаусса использовать выбор опорного элемента: об этом можно почитать в источнике кода @@ -171,7 +172,6 @@ namespace { // * (простое описание для понимания) // * [3] http://ikrisoft.blogspot.com/2015/01/ransac-with-contrario-approach.html - constexpr float min_dst2 = 1e-6f; const int n_matches = points_lhs.size(); // https://en.wikipedia.org/wiki/Random_sample_consensus#Parameters @@ -189,10 +189,10 @@ namespace { randomSample(sample, n_matches, n_samples, &seed, [&sample, &points_lhs, &points_rhs](int x) { for (int prev : sample) { auto left_vec = (points_lhs[prev] - points_lhs[x]); - if (std::hypot(left_vec.x, left_vec.y) < min_dst2) + if (std::hypot(left_vec.x, left_vec.y) < MIN_DST2) return false; auto right_vec = (points_rhs[prev] - points_rhs[x]); - if (std::hypot(right_vec.x, right_vec.y) < min_dst2) + if (std::hypot(right_vec.x, right_vec.y) < MIN_DST2) return false; } return true;