From 4e7eea5ac21228f917cd5e486faeaed71fa354c2 Mon Sep 17 00:00:00 2001 From: Egor Date: Fri, 6 Mar 2026 22:04:12 +0300 Subject: [PATCH 1/6] done --- src/phg/sift/sift.cpp | 165 +++++++++++++++++++++++++----------------- tests/test_sift.cpp | 2 +- 2 files changed, 100 insertions(+), 67 deletions(-) diff --git a/src/phg/sift/sift.cpp b/src/phg/sift/sift.cpp index 72047711..297dd90d 100755 --- a/src/phg/sift/sift.cpp +++ b/src/phg/sift/sift.cpp @@ -107,17 +107,19 @@ std::vector phg::buildOctaves(const cv::Mat& img, const phg:: // для простоты в каждой октаве будем каждый раз блюрить базовую картинку с полной сигмой // можно подумать, как сделать эффективнее - для построения n+1 слоя доблюревать уже поблюренный n-ый слой, так чтобы в итоге получилась такая же сигма // это будет немного быстрее, тк нужно более маленькое ядро свертки на каждый шаг + double k = std::pow(2., 1. / s); 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_p = std::pow(k, (double)(i - 1)) * sigma0; + double sigma_layer = sigma_p * k; + // вычтем sigma0 чтобы размыть ровно до нужной суммарной сигмы + sigma_layer = std::sqrt(sigma_layer * sigma_layer - sigma_p * sigma_p); + cv::GaussianBlur(oct.layers[i - 1], oct.layers[i], cv::Size(), sigma_layer, sigma_layer); } // подготавливаем базовый слой для следующей октавы if (o + 1 < n_octaves) { // используется в opencv, формула для пересчета ключевых точек: pt_upscaled = 2^o * pt_downscaled - // TODO cv::resize(даунскейлим текущий слой в два раза, без интерполяции, просто сабсепмлинг); + base = downsample2x(oct.layers[s]); // можно использовать и downsample2x_avg(oct.layers[s]), это позволяет потом заапскейлить слои обратно до оригинального разрешения без сдвига // но потребуется везде изменить формулу для пересчета ключевых точек: pt_upscaled = (pt_downscaled + 0.5) * 2^o - 0.5 @@ -137,7 +139,9 @@ std::vector phg::buildDoG(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 = nL.at(yi, xi) + pL.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) - pL.at(yi + 1, xi) - nL.at(yi - 1, xi) + pL.at(yi - 1, xi)) * 0.25f; cv::Matx33f H(dxx, dxy, dxs, dxy, dyy, dys, dxs, dys, dss); @@ -273,10 +302,10 @@ std::vector phg::findScaleSpaceExtrema(const std::vector phg::findScaleSpaceExtrema(const std::vector (r + 1) * (r + 1) / r) + break; } // скейлим координаты точек обратно до родных размеров картинки @@ -379,39 +408,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 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 = TODO; -// float angle = std::atan2(TODO); // [-pi, pi] + 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 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 weight = std::exp(-(dx * dx + dy * dy) / (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 bin = angle_deg / 360 * 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; -// } + float frac = bin - bin0; + if (!params.enable_orientation_bin_interpolation) { + frac = 0.f; + } // -// histogram[bin0] += TODO; -// histogram[bin1] += TODO; + histogram[bin0] += (1.f - frac) * mag * weight; + histogram[bin1] += frac * weight * mag; } } @@ -450,20 +479,20 @@ 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 offset = 0.5 * (left - right) / (left + right - 2. * center); + 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 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; + float angle = bin_real * 360.f / n_bins; // -// cv::KeyPoint new_kp = kp; -// new_kp.angle = angle; -// oriented_kpts.push_back(new_kp); + cv::KeyPoint new_kp = kp; + new_kp.angle = angle; + oriented_kpts.push_back(new_kp); } } } @@ -529,7 +558,7 @@ std::pair> phg::computeDescriptors(const std: float sin_a = std::sin(kp_angle_rad); // для гауссового взвешивания: затухающий вклад градиентов с краев картинки - float sigma_desc = (float)n_spatial_bins * 0.5f; + float sigma_desc = (float)n_spatial_bins * 0.4f; std::vector desc(n_dims, 0.f); @@ -570,15 +599,17 @@ std::pair> phg::computeDescriptors(const std: // подсчет бина в гистограммке градиентов внутри пространственного бина float bin_o = angle_invariant * n_orient_bins / CV_2PI; + if (bin_o < 0.f) + bin_o += n_orient_bins; if (bin_o >= n_orient_bins) 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 +640,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] += wy * wx * wo * weighted_mag; } } } @@ -621,8 +652,8 @@ std::pair> phg::computeDescriptors(const std: 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; } } } @@ -783,3 +814,5 @@ void phg::SIFT::savePyramid(const std::string& name, const std::vector& } } } + + diff --git a/tests/test_sift.cpp b/tests/test_sift.cpp index cf3bd7df..7433f9bb 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 From 32e6e57c097339cb817cf1b2af211e8f3e84b94e Mon Sep 17 00:00:00 2001 From: EgorM2718 Date: Thu, 19 Mar 2026 21:01:20 +0300 Subject: [PATCH 2/6] main part done --- src/phg/matching/descriptor_matcher.cpp | 98 ++++++++++++++++-------- src/phg/matching/flann_matcher.cpp | 26 ++++++- src/phg/sfm/homography.cpp | 99 ++++++++++++++----------- src/phg/sfm/panorama_stitcher.cpp | 25 ++++++- 4 files changed, 168 insertions(+), 80 deletions(-) diff --git a/src/phg/matching/descriptor_matcher.cpp b/src/phg/matching/descriptor_matcher.cpp index f4bcd871..92d59f11 100644 --- a/src/phg/matching/descriptor_matcher.cpp +++ b/src/phg/matching/descriptor_matcher.cpp @@ -7,8 +7,21 @@ void phg::DescriptorMatcher::filterMatchesRatioTest(const std::vector &filtered_matches) { filtered_matches.clear(); - - throw std::runtime_error("not implemented yet"); + + const float ratio_threshold = 0.75f; + + for (const auto& match_pair : matches) { + if (match_pair.size() < 2) continue; + + float dist1 = match_pair[0].distance; + float dist2 = match_pair[1].distance; + + if (dist2 == 0) continue; + + if (dist1 < ratio_threshold * dist2) { + filtered_matches.push_back(match_pair[0]); + } + } } @@ -32,45 +45,70 @@ void phg::DescriptorMatcher::filterMatchesClusters(const std::vector cv::Mat points_query(n_matches, 2, CV_32FC1); cv::Mat points_train(n_matches, 2, CV_32FC1); for (int i = 0; i < n_matches; ++i) { - points_query.at(i) = keypoints_query[matches[i].queryIdx].pt; - points_train.at(i) = keypoints_train[matches[i].trainIdx].pt; + points_query.at(i, 0) = keypoints_query[matches[i].queryIdx].pt.x; + points_query.at(i, 1) = keypoints_query[matches[i].queryIdx].pt.y; + points_train.at(i, 0) = keypoints_train[matches[i].trainIdx].pt.x; + points_train.at(i, 1) = keypoints_train[matches[i].trainIdx].pt.y; } // // // размерность всего 2, так что точное KD-дерево -// std::shared_ptr index_params = flannKdTreeIndexParams(TODO); -// std::shared_ptr search_params = flannKsTreeSearchParams(TODO); + std::shared_ptr index_params = flannKdTreeIndexParams(2); + 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); + 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); + 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); + 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; -// } + 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 + for (int i = 0; i < n_matches; ++i) { + std::vector valid_query_neighbors; + for (int j = 0; j < total_neighbours; ++j) { + if (distances2_query.at(i, j) <= radius2_query) { + valid_query_neighbors.push_back(indices_query.at(i, j)); + } + } + + int shared_count = 0; + for (int j = 0; j < total_neighbours; ++j) { + if (distances2_train.at(i, j) <= radius2_train) { + int train_neighbor_idx = indices_train.at(i, j); + + if (std::find(valid_query_neighbors.begin(), valid_query_neighbors.end(), train_neighbor_idx) != valid_query_neighbors.end()) { + shared_count++; + } + } + } + + if (shared_count >= consistent_matches) { + filtered_matches.push_back(matches[i]); + } +} } diff --git a/src/phg/matching/flann_matcher.cpp b/src/phg/matching/flann_matcher.cpp index 9e9f5180..b87ad357 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,25 @@ 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"); + matches.clear(); + + if (query_desc.empty() || !flann_index) { + return; + } + + cv::Mat indices, distances; + flann_index->knnSearch(query_desc, indices, distances, 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 trainIdx = indices.at(i, j); + float dist = std::sqrt(distances.at(i, j)); + + if (trainIdx >= 0) { + matches[i].push_back(cv::DMatch(i, trainIdx, dist)); + } + } + } } diff --git a/src/phg/sfm/homography.cpp b/src/phg/sfm/homography.cpp index 5cbc780c..2e67b189 100644 --- a/src/phg/sfm/homography.cpp +++ b/src/phg/sfm/homography.cpp @@ -84,8 +84,9 @@ 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, -1, 0, 0, 0, x0*x1, y0*x1, -x1}); + A.push_back({0, 0, 0, -x0, -y0, -1, x0*y1, y0*y1, -y1}); + } int res = gauss(A, H); @@ -168,57 +169,57 @@ namespace { // * (простое описание для понимания) // * [3] http://ikrisoft.blogspot.com/2015/01/ransac-with-contrario-approach.html -// const int n_matches = points_lhs.size(); + const int n_matches = points_lhs.size(); // // // https://en.wikipedia.org/wiki/Random_sample_consensus#Parameters -// const int n_trials = TODO; + const int n_trials = 5000; // -// const int n_samples = TODO; -// uint64_t seed = 1; -// const double reprojection_error_threshold_px = 2; + const int n_samples = 4; + uint64_t seed = 1; + const double reprojection_error_threshold_px = 2; // -// int best_support = 0; -// cv::Mat best_H; + 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::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; + } + } // -// std::cout << "estimateHomographyRANSAC : support: " << best_support << "/" << n_matches << std::endl; + if (support > best_support) { + best_support = support; + best_H = H; // -// if (best_support == n_matches) { -// break; -// } -// } -// } + std::cout << "estimateHomographyRANSAC : support: " << best_support << "/" << n_matches << std::endl; // -// std::cout << "estimateHomographyRANSAC : best support: " << best_support << "/" << n_matches << std::endl; + if (best_support == n_matches) { + break; + } + } + } // -// if (best_support == 0) { -// throw std::runtime_error("estimateHomographyRANSAC : failed to estimate homography"); -// } +// std::cout << "estimateHomographyRANSAC : best support: " << best_support << "/" << n_matches << std::endl; // -// return best_H; + if (best_support == 0) { + throw std::runtime_error("estimateHomographyRANSAC : failed to estimate homography"); + } + + return best_H; } } @@ -238,7 +239,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"); + double x = pt.x, y = pt.y; + double w = T.at(2,0) * x + T.at(2,1) * y + T.at(2,2); + + if (std::abs(w) < 1e-8) { + throw std::runtime_error("transformPoint: division by zero"); + } + + return cv::Point2d( + (T.at(0,0) * x + T.at(0,1) * y + T.at(0,2)) / w, + (T.at(1,0) * x + T.at(1,1) * y + T.at(1,2)) / 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..c5f4101d 100644 --- a/src/phg/sfm/panorama_stitcher.cpp +++ b/src/phg/sfm/panorama_stitcher.cpp @@ -21,9 +21,28 @@ cv::Mat phg::stitchPanorama(const std::vector &imgs, // вектор гомографий, для каждой картинки описывает преобразование до корня std::vector Hs(n_images); { - // здесь надо посчитать вектор Hs - // при этом можно обойтись n_images - 1 вызовами функтора homography_builder - throw std::runtime_error("not implemented yet"); + std::vector computed(n_images, false); + std::function compute_H = [&](int i) { + if (computed[i]) return; + + if (parent[i] == -1) { + Hs[i] = cv::Mat::eye(3, 3, CV_64F); + computed[i] = true; + return; + } + + compute_H(parent[i]); + + cv::Mat H_local = homography_builder(imgs[i], imgs[parent[i]]); + + Hs[i] = Hs[parent[i]] * H_local; + + computed[i] = true; + }; + + for (int i = 0; i < n_images; ++i) { + compute_H(i); + } } bbox2 bbox; From f63a41a759cfc024b802e68355f0cb51bd5ee8f4 Mon Sep 17 00:00:00 2001 From: EgorM2718 Date: Thu, 19 Mar 2026 21:50:25 +0300 Subject: [PATCH 3/6] a contrario done --- src/phg/matching/descriptor_matcher.cpp | 1 - src/phg/sfm/homography.cpp | 143 ++++++++++++++---------- 2 files changed, 87 insertions(+), 57 deletions(-) diff --git a/src/phg/matching/descriptor_matcher.cpp b/src/phg/matching/descriptor_matcher.cpp index 92d59f11..bc30a202 100644 --- a/src/phg/matching/descriptor_matcher.cpp +++ b/src/phg/matching/descriptor_matcher.cpp @@ -87,7 +87,6 @@ void phg::DescriptorMatcher::filterMatchesClusters(const std::vector } // // метч остается, если левое и правое множества первых total_neighbors соседей в радиусах поиска(radius2_query, radius2_train) имеют как минимум consistent_matches общих элементов -// // TODO заполнить filtered_matches for (int i = 0; i < n_matches; ++i) { std::vector valid_query_neighbors; for (int j = 0; j < total_neighbours; ++j) { diff --git a/src/phg/sfm/homography.cpp b/src/phg/sfm/homography.cpp index 2e67b189..4b57bb84 100644 --- a/src/phg/sfm/homography.cpp +++ b/src/phg/sfm/homography.cpp @@ -1,5 +1,5 @@ #include "homography.h" - +#include #include #include @@ -158,70 +158,101 @@ namespace { } cv::Mat estimateHomographyRANSAC(const std::vector &points_lhs, const std::vector &points_rhs) - { - if (points_lhs.size() != points_rhs.size()) { - throw std::runtime_error("findHomography: points_lhs.size() != points_rhs.size()"); +{ + if (points_lhs.size() != points_rhs.size()) { + throw std::runtime_error("findHomography: points_lhs.size() != points_rhs.size()"); + } + + const int n_matches = points_lhs.size(); + if (n_matches < 4) { + throw std::runtime_error("estimateHomographyRANSAC: not enough points"); + } + + if (n_matches == 4) { + return estimateHomography4Points(points_lhs[0], points_lhs[1], points_lhs[2], points_lhs[3], + points_rhs[0], points_rhs[1], points_rhs[2], points_rhs[3]); + } + + double min_x = points_rhs[0].x, max_x = points_rhs[0].x; + double min_y = points_rhs[0].y, max_y = points_rhs[0].y; + for (const auto& p : points_rhs) { + min_x = std::min(min_x, (double)p.x); max_x = std::max(max_x, (double)p.x); + min_y = std::min(min_y, (double)p.y); max_y = std::max(max_y, (double)p.y); + } + double target_area = (max_x - min_x) * (max_y - min_y); + if (target_area < 1.0) target_area = 1e6; + + const int n_trials = 5000; + const int n_samples = 4; + uint64_t seed = 1; + + double best_log_nfa = std::numeric_limits::max(); + cv::Mat best_H; + int best_support = 0; + + double log_N_minus_4 = std::log(std::max(1, n_matches - 4)); + double lgamma_n_plus_1 = std::lgamma(n_matches + 1.0); + double lgamma_4_plus_1 = std::lgamma(n_samples + 1.0); + + std::vector sample; + std::vector errors(n_matches); + + for (int i_trial = 0; i_trial < n_trials; ++i_trial) { + randomSample(sample, n_matches, n_samples, &seed); + + cv::Mat H; + try { + 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]]); + } catch (...) { + continue; } - // TODO Дополнительный балл, если вместо обычной версии будет использована модификация a-contrario RANSAC - // * [1] Automatic Homographic Registration of a Pair of Images, with A Contrario Elimination of Outliers. (Lionel Moisan, Pierre Moulon, Pascal Monasse) - // * [2] Adaptive Structure from Motion with a contrario model estimation. (Pierre Moulon, Pascal Monasse, Renaud Marlet) - // * (простое описание для понимания) - // * [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 = 5000; -// - 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); - - 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; - } + for (int i = 0; i < n_matches; ++i) { + try { + cv::Point2d proj = phg::transformPoint(points_lhs[i], H); + cv::Point2d diff = proj - cv::Point2d(points_rhs[i]); + errors[i] = diff.x * diff.x + diff.y * diff.y; + } catch (...) { + errors[i] = std::numeric_limits::max(); } -// - if (support > best_support) { - best_support = support; + } + + std::vector sorted_errors = errors; + std::sort(sorted_errors.begin(), sorted_errors.end()); + + for (int k = 5; k <= n_matches; ++k) { + double e2 = sorted_errors[k - 1]; + + if (e2 >= std::numeric_limits::max() / 2.0) break; + + double alpha = (CV_PI * e2) / target_area; + if (alpha >= 1.0) continue; + alpha = std::max(alpha, 1e-12); + + double lgamma_k_plus_1 = std::lgamma(k + 1.0); + double log_binom_n_k = lgamma_n_plus_1 - lgamma_k_plus_1 - std::lgamma(n_matches - k + 1.0); + double log_binom_k_4 = lgamma_k_plus_1 - lgamma_4_plus_1 - std::lgamma(k - n_samples + 1.0); + + double log_nfa = log_N_minus_4 + log_binom_n_k + log_binom_k_4 + (k - n_samples) * std::log(alpha); + + if (log_nfa < best_log_nfa) { + best_log_nfa = log_nfa; best_H = H; -// - std::cout << "estimateHomographyRANSAC : support: " << best_support << "/" << n_matches << std::endl; -// - if (best_support == n_matches) { - break; - } + best_support = k; } } -// -// 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; + if (best_log_nfa > 0) { + throw std::runtime_error("estimateHomographyRANSAC : failed to estimate homography"); } + std::cout << "estimateHomographyRANSAC : support: " << best_support << "/" << n_matches << std::endl; + + return best_H; +} + } cv::Mat phg::findHomography(const std::vector &points_lhs, const std::vector &points_rhs) From d356bd1b7d745335189a47dae16c2d30e5fb9293 Mon Sep 17 00:00:00 2001 From: EgorM2718 Date: Thu, 19 Mar 2026 21:51:41 +0300 Subject: [PATCH 4/6] done --- tests/test_matching.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_matching.cpp b/tests/test_matching.cpp index adaac65e..c9e927e1 100644 --- a/tests/test_matching.cpp +++ b/tests/test_matching.cpp @@ -20,7 +20,7 @@ // TODO enable both toggles for testing custom detector & matcher #define ENABLE_MY_DESCRIPTOR 0 -#define ENABLE_MY_MATCHING 0 +#define ENABLE_MY_MATCHING 1 #define ENABLE_GPU_BRUTEFORCE_MATCHER 0 // TODO disable for local testing but do not commit From 038db8def59adfec9ac24cb34ba18616a5852a90 Mon Sep 17 00:00:00 2001 From: EgorM2718 Date: Fri, 20 Mar 2026 20:49:40 +0300 Subject: [PATCH 5/6] - --- src/phg/matching/flann_matcher.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/phg/matching/flann_matcher.cpp b/src/phg/matching/flann_matcher.cpp index b87ad357..114ccfe2 100644 --- a/src/phg/matching/flann_matcher.cpp +++ b/src/phg/matching/flann_matcher.cpp @@ -7,7 +7,7 @@ phg::FlannMatcher::FlannMatcher() { // параметры для приближенного поиска index_params = flannKdTreeIndexParams(4); - search_params = flannKsTreeSearchParams(32); + search_params = flannKsTreeSearchParams(16); } void phg::FlannMatcher::train(const cv::Mat &train_desc) From 70d536531c2618780c09abd495036d2ffd9839b5 Mon Sep 17 00:00:00 2001 From: EgorM2718 Date: Fri, 20 Mar 2026 20:59:41 +0300 Subject: [PATCH 6/6] -- --- src/phg/matching/flann_matcher.cpp | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/src/phg/matching/flann_matcher.cpp b/src/phg/matching/flann_matcher.cpp index 114ccfe2..9c20857d 100644 --- a/src/phg/matching/flann_matcher.cpp +++ b/src/phg/matching/flann_matcher.cpp @@ -7,7 +7,7 @@ phg::FlannMatcher::FlannMatcher() { // параметры для приближенного поиска index_params = flannKdTreeIndexParams(4); - search_params = flannKsTreeSearchParams(16); + search_params = flannKsTreeSearchParams(32); } void phg::FlannMatcher::train(const cv::Mat &train_desc) @@ -23,18 +23,23 @@ void phg::FlannMatcher::knnMatch(const cv::Mat &query_desc, std::vectorknnSearch(query_desc, indices, distances, k, *search_params); matches.resize(query_desc.rows); for (int i = 0; i < query_desc.rows; ++i) { + matches[i].reserve(k); + + const int* indices_ptr = indices.ptr(i); + const float* distances_ptr = distances.ptr(i); + for (int j = 0; j < k; ++j) { - int trainIdx = indices.at(i, j); - float dist = std::sqrt(distances.at(i, j)); - + int trainIdx = indices_ptr[j]; if (trainIdx >= 0) { - matches[i].push_back(cv::DMatch(i, trainIdx, dist)); + matches[i].emplace_back(i, trainIdx, std::sqrt(distances_ptr[j])); } } }