diff --git a/.gitignore b/.gitignore index 54b965f..dac4730 100755 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,5 @@ .idea build cmake-build* +.cache +data/debug \ No newline at end of file diff --git a/src/phg/matching/cl/bruteforce_matcher.cl b/src/phg/matching/cl/bruteforce_matcher.cl index e51a34a..4e77b11 100644 --- a/src/phg/matching/cl/bruteforce_matcher.cl +++ b/src/phg/matching/cl/bruteforce_matcher.cl @@ -22,6 +22,7 @@ __kernel void bruteforce_matcher(__global const float* train, // храним KEYPOINTS_PER_WG=4 дескриптора-query: __local float query_local[KEYPOINTS_PER_WG * NDIM]; + __local float dist2_for_reduction[NDIM]; // храним два лучших сопоставления для каждого дескриптора-query: __local uint res_train_idx_local[KEYPOINTS_PER_WG * 2]; __local float res_distance2_local[KEYPOINTS_PER_WG * 2]; // храним квадраты чтобы не считать корень до самого последнего момента @@ -33,6 +34,11 @@ __kernel void bruteforce_matcher(__global const float* train, // грузим 4 дескриптора-query (для каждого из четырех дескрипторов каждый поток грузит значение своей размерности dim_id) // TODO: т.е. надо прогрузить в query_local все KEYPOINTS_PER_WG=4 дескриптора из query (начиная с индекса query_id0) (а если часть из них выходит за пределы n_query_desc - грузить нули) + for (int query_local_i = 0; query_local_i < KEYPOINTS_PER_WG; query_local_i++) { + int query_id = query_id0 + query_local_i; + query_local[query_local_i * NDIM + dim_id] = query_id < n_query_desc ? query[query_id * NDIM + dim_id] : 0.f; + } + barrier(CLK_LOCAL_MEM_FENCE); // дожидаемся прогрузки наших дескрипторов-запросов for (int train_idx = 0; train_idx < n_train_desc; ++train_idx) { @@ -43,13 +49,15 @@ __kernel void bruteforce_matcher(__global const float* train, // до дескриптора-train в глобальной памяти (#train_idx) // TODO посчитать квадрат расстояния по нашей размерности (dim_id) и сохранить его в нашу ячейку в dist2_for_reduction + float d = query_local[query_local_i * NDIM + dim_id] - train_value_dim; + dist2_for_reduction[dim_id] = d * d; barrier(CLK_LOCAL_MEM_FENCE); // TODO суммируем редукцией все что есть в dist2_for_reduction int step = NDIM / 2; while (step > 0) { if (dim_id < step) { - // TODO + dist2_for_reduction[dim_id] += dist2_for_reduction[dim_id + step]; } barrier(CLK_LOCAL_MEM_FENCE); step /= 2; @@ -63,18 +71,24 @@ __kernel void bruteforce_matcher(__global const float* train, #define SECOND_BEST_INDEX 1 // пытаемся улучшить самое лучшее сопоставление для локального дескриптора - if (dist2 <= res_distance2_local[query_local_i * 2 + BEST_INDEX]) { + if (dist2 < res_distance2_local[query_local_i * 2 + BEST_INDEX]) { // не забываем что прошлое лучшее сопоставление теперь стало вторым по лучшевизне (на данный момент) res_distance2_local[query_local_i * 2 + SECOND_BEST_INDEX] = res_distance2_local[query_local_i * 2 + BEST_INDEX]; res_train_idx_local[query_local_i * 2 + SECOND_BEST_INDEX] = res_train_idx_local[query_local_i * 2 + BEST_INDEX]; // TODO заменяем нашим (dist2, train_idx) самое лучшее сопоставление для локального дескриптора - } else if (dist2 <= res_distance2_local[query_local_i * 2 + SECOND_BEST_INDEX]) { // может мы улучшили хотя бы второе по лучшевизне сопоставление? + res_distance2_local[query_local_i * 2 + BEST_INDEX] = dist2; + res_train_idx_local[query_local_i * 2 + BEST_INDEX] = train_idx; + } else if (dist2 < res_distance2_local[query_local_i * 2 + SECOND_BEST_INDEX]) { // может мы улучшили хотя бы второе по лучшевизне сопоставление? // TODO заменяем второе по лучшевизне сопоставление для локального дескриптора + res_distance2_local[query_local_i * 2 + SECOND_BEST_INDEX] = dist2; + res_train_idx_local[query_local_i * 2 + SECOND_BEST_INDEX] = train_idx; } } } } + barrier(CLK_LOCAL_MEM_FENCE); + // итак, мы нашли два лучших сопоставления для наших KEYPOINTS_PER_WG дескрипторов, надо сохрнить эти результаты в глобальную память if (dim_id < KEYPOINTS_PER_WG * 2) { // полагаемся на то что нам надо прогрузить KEYPOINTS_PER_WG*2==4*2 #include "flann_factory.h" +#include void phg::DescriptorMatcher::filterMatchesRatioTest(const std::vector> &matches, std::vector &filtered_matches) { filtered_matches.clear(); + const float threshold2 = 0.705f * 0.705f; - throw std::runtime_error("not implemented yet"); + for (const auto& knn : matches) { + if (knn.size() >= 2 && knn[0].distance < threshold2 * knn[1].distance) { + filtered_matches.push_back(knn[0]); + } + } } @@ -35,42 +41,65 @@ 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; -// } -// + + // размерность всего 2, так что точное KD-дерево + std::shared_ptr index_params = flannKdTreeIndexParams(1); + 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_neighbors соседей в радиусах поиска(radius2_query, radius2_train) имеют как минимум consistent_matches общих элементов // // TODO заполнить filtered_matches + + std::unordered_set query_neighbors; + for (size_t matchIdx = 0; matchIdx < n_matches; matchIdx++) { + for (size_t neighborIdx = 0; neighborIdx < total_neighbours; neighborIdx++) { + if (distances2_query.at(matchIdx, neighborIdx) <= radius2_query) { + query_neighbors.insert(indices_query.at(matchIdx, neighborIdx)); + } + } + + size_t count = 0; + for (size_t neighborIdx = 0; neighborIdx < total_neighbours; neighborIdx++) { + if (distances2_train.at(matchIdx, neighborIdx) <= radius2_train && + query_neighbors.find(indices_train.at(matchIdx, neighborIdx)) != query_neighbors.end()) { + count++; + } + } + + if (count >= consistent_matches) { + filtered_matches.emplace_back(matches[matchIdx]); + } + + query_neighbors.clear(); + } } diff --git a/src/phg/matching/flann_matcher.cpp b/src/phg/matching/flann_matcher.cpp index 9e9f518..a0267a2 100644 --- a/src/phg/matching/flann_matcher.cpp +++ b/src/phg/matching/flann_matcher.cpp @@ -1,13 +1,15 @@ #include #include "flann_matcher.h" #include "flann_factory.h" +#include "opencv2/core/hal/interface.h" +#include "opencv2/core/types.hpp" 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 +19,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"); + cv::Mat indices(query_desc.rows, k, CV_32SC1); + cv::Mat distances2(query_desc.rows, k, CV_32FC1); + + flann_index->knnSearch( + query_desc, + indices, + distances2, + k, + *search_params + ); + + matches.resize(query_desc.rows); + for (size_t rowIdx = 0; rowIdx < query_desc.rows; rowIdx++) { + for (size_t idx = 0; idx < k; idx++) { + matches[rowIdx].emplace_back( + rowIdx, + indices.at(rowIdx, idx), + distances2.at(rowIdx, idx) + ); + } + } } diff --git a/src/phg/sfm/homography.cpp b/src/phg/sfm/homography.cpp index 5cbc780..2dd6107 100644 --- a/src/phg/sfm/homography.cpp +++ b/src/phg/sfm/homography.cpp @@ -2,6 +2,7 @@ #include #include +#include namespace { @@ -84,8 +85,8 @@ namespace { double w1 = ws1[i]; // 8 elements of matrix + free term as needed by gauss routine -// A.push_back({TODO}); -// A.push_back({TODO}); + A.push_back({0.0, 0.0, 0.0, -x0 * w1, -y0 * w1, -w0 * w1, x0 * y1, y0 * y1, -w0 * y1}); + A.push_back({x0 * w1, y0 * w1, w0 * w1, 0.0, 0.0, 0.0, -x0 * x1, -y0 * x1, w0 * x1}); } int res = gauss(A, H); @@ -168,57 +169,59 @@ 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 double p = 0.99; + const double w = 0.5; + const int n_trials = std::round(std::log(1.0 - p) / std::log(1 - std::pow(w, 4.0))); + + 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; + } + } + + 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 +241,11 @@ 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"); + cv::Mat out = T * cv::Mat(cv::Point3d(pt.x, pt.y, 1.0)); + return { + out.at(0) / out.at(2), + out.at(1) / out.at(2) + }; } 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..b80ed51 100644 --- a/src/phg/sfm/panorama_stitcher.cpp +++ b/src/phg/sfm/panorama_stitcher.cpp @@ -4,6 +4,34 @@ #include #include +cv::Mat buildTransform( + size_t imgIdx, + std::vector& used, + std::vector& Hs, + const std::vector &imgs, + const std::vector &parent, + std::function &homography_builder +) { + if (used[imgIdx]) { + return Hs[imgIdx]; + } + + used[imgIdx] = true; + + if (parent[imgIdx] == -1) { + return Hs[imgIdx] = cv::Mat::eye(3, 3, CV_64FC1); + } + + return Hs[imgIdx] = homography_builder(imgs[imgIdx], imgs[parent[imgIdx]]) * buildTransform( + parent[imgIdx], + used, + Hs, + imgs, + parent, + homography_builder + ); +} + /* * imgs - список картинок * parent - список индексов, каждый индекс указывает, к какой картинке должна быть приклеена текущая картинка @@ -23,7 +51,17 @@ cv::Mat phg::stitchPanorama(const std::vector &imgs, { // здесь надо посчитать вектор Hs // при этом можно обойтись n_images - 1 вызовами функтора homography_builder - throw std::runtime_error("not implemented yet"); + std::vector used(n_images, false); + for (size_t imgIdx = 0; imgIdx < n_images; ++imgIdx) { + buildTransform( + imgIdx, + used, + Hs, + imgs, + parent, + homography_builder + ); + } } bbox2 bbox; @@ -46,19 +84,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 = transformPoint(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(); }); diff --git a/src/phg/sift/sift.cpp b/src/phg/sift/sift.cpp index 7204771..9f5fcb9 100755 --- a/src/phg/sift/sift.cpp +++ b/src/phg/sift/sift.cpp @@ -1,8 +1,12 @@ #include "sift.h" #include "libutils/rasserts.h" +#include +#include #include #include +#include +#include #include #include #include @@ -79,6 +83,7 @@ cv::Mat phg::toGray32F(const cv::Mat& img) std::vector phg::buildOctaves(const cv::Mat& img, const phg::SIFTParams& p, int verbose_level) { const int s = p.n_octave_layers; + const double sth_root_two = std::pow(2.0, 1.0 / s); const double sigma0 = p.sigma; // взятое с потолка значение блюра который уже есть в картинке. используем для того, чтобы не так сильно блюрить базовую картинку и не терять лишний раз фичи // upd: хотя llm не соглашается со "взятое с потолка": @@ -107,18 +112,23 @@ std::vector phg::buildOctaves(const cv::Mat& img, const phg:: // для простоты в каждой октаве будем каждый раз блюрить базовую картинку с полной сигмой // можно подумать, как сделать эффективнее - для построения n+1 слоя доблюревать уже поблюренный n-ый слой, так чтобы в итоге получилась такая же сигма // это будет немного быстрее, тк нужно более маленькое ядро свертки на каждый шаг + double factor = sth_root_two; 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 * factor; + sigma_layer = std::sqrt(sigma_layer * sigma_layer - sigma0 * sigma0); + cv::GaussianBlur(oct.layers[0], oct.layers[i], cv::Size(), sigma_layer, sigma_layer); + factor *= sth_root_two; } // подготавливаем базовый слой для следующей октавы if (o + 1 < n_octaves) { // используется в opencv, формула для пересчета ключевых точек: pt_upscaled = 2^o * pt_downscaled // TODO 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 @@ -139,6 +149,9 @@ std::vector phg::buildDoG(const std::vector phg::findScaleSpaceExtrema(const std::vector curr_offsets = {-1, 0, 1}) { + for (const float* layer : {prev, next}) { + for (int offset : {-1, 0, 1}) { + check(layer[x + offset]); + } + } - // TODO проверить локальный максимум на текущем скейле + for (int offset : curr_offsets) { + check(curr[x + offset]); + } + }; + // TODO проверить локальный максимум на текущем скейле + checkNeighbors(cp, c, cn, {-1, 1}); if (!is_max && !is_min) continue; // TODO проверить локальный максимум на предыдущем скейле - + checkNeighbors(pp, p, pn); if (!is_max && !is_min) continue; // TODO проверить локальный максимум на следующем скейле - + checkNeighbors(np, n, nn); if (!is_max && !is_min) continue; @@ -238,13 +263,12 @@ 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) - 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); @@ -288,6 +312,17 @@ std::vector phg::findScaleSpaceExtrema(const std::vector (r + 1) * (r + 1) * det) { + break; + } } // скейлим координаты точек обратно до родных размеров картинки @@ -379,39 +414,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::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(-(dx * dx + dy * dy) / (2.f * sigma_win * sigma_win)); + if (!params.enable_orientation_gaussian_weighting) { + weight = 1.f; + } + + // голосуем в гистограмме направлений. находим два ближайших бина и гладко распределяем голос между ними + // в таком случае, голос попавший близко к границе между бинами, проголосует поровну за оба бина + float bin = angle_deg * n_bins / 360.f; + 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] += weight * mag * (1.f - frac); + histogram[bin1] += weight * mag * frac; } } @@ -441,29 +476,29 @@ 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); + // хотим найти угол дескриптора точнее = зафитить параболу по трем точкам (i-1, left), (i, center), (i+1, right) + // у параболы f(x) = ax^2 + bx + c, экстремум в точке x = offset = -b/(2a) + // f(-1) = left, f(0) = center, f(1) = right + // f(0) = c = center + // f(1) = a + b + c = right + // f(-1) = a - b + c = left + // f(1) + f(-1) = 2a + 2c -> a = (left + right - 2 * center) / 2 + // f(1) - f(-1) = 2b -> b = (right - left) / 2 + + float offset = -0.5f * (right - left) / (left + right - 2.f * 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 angle = bin_real * 360.f / n_bins; + + cv::KeyPoint new_kp = kp; + new_kp.angle = angle; + oriented_kpts.push_back(new_kp); } } } @@ -574,11 +609,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 +644,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; } } } @@ -621,8 +656,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; } } } diff --git a/tests/test_matching.cpp b/tests/test_matching.cpp index adaac65..26889de 100644 --- a/tests/test_matching.cpp +++ b/tests/test_matching.cpp @@ -1,3 +1,4 @@ +#include #include #include @@ -19,9 +20,9 @@ // TODO enable both toggles for testing custom detector & matcher -#define ENABLE_MY_DESCRIPTOR 0 -#define ENABLE_MY_MATCHING 0 -#define ENABLE_GPU_BRUTEFORCE_MATCHER 0 +#define ENABLE_MY_DESCRIPTOR 1 +#define ENABLE_MY_MATCHING 1 +#define ENABLE_GPU_BRUTEFORCE_MATCHER 1 // TODO disable for local testing but do not commit #define SERVER_TESTING 1 diff --git a/tests/test_sift.cpp b/tests/test_sift.cpp index cf3bd7d..7433f9b 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