diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml index 66e2997..c1c9caf 100755 --- a/.github/workflows/cmake.yml +++ b/.github/workflows/cmake.yml @@ -2,6 +2,7 @@ name: CMake on: [push, pull_request] + env: BUILD_TYPE: RelWithDebInfo diff --git a/.gitignore b/.gitignore index 54b965f..e033a90 100755 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,5 @@ .idea build cmake-build* +*.zip +opencv-* \ No newline at end of file diff --git a/src/phg/matching/descriptor_matcher.cpp b/src/phg/matching/descriptor_matcher.cpp index f4bcd87..e8c1060 100644 --- a/src/phg/matching/descriptor_matcher.cpp +++ b/src/phg/matching/descriptor_matcher.cpp @@ -3,25 +3,28 @@ #include #include "flann_factory.h" -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"); + filtered_matches.reserve(matches.size()); + for (const auto& match: matches) { + if (match.size() < 2) continue; + if (match[0].distance < 0.75f * match[1].distance) + filtered_matches.push_back(match[0]); + } } -void phg::DescriptorMatcher::filterMatchesClusters(const std::vector &matches, +void phg::DescriptorMatcher::filterMatchesClusters(const std::vector& matches, const std::vector keypoints_query, const std::vector keypoints_train, - std::vector &filtered_matches) -{ + 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 +38,51 @@ 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 + filtered_matches.reserve(n_matches); + for (int i = 0; i < n_matches; ++i) { + std::set s; + for (int j = 0; j < total_neighbours && distances2_query.at(i, j) <= radius2_query; ++j) + s.insert(indices_query.at(i, j)); + int intersections = 0; + for (int j = 0; j < total_neighbours && distances2_train.at(i, j) <= radius2_train; ++j) + if (s.count(indices_train.at(i, j))) ++intersections; + if (intersections >= 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 9e9f518..c3484f5 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,13 @@ 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, dst; + flann_index->knnSearch(query_desc, idx, dst, k, *search_params); + matches.resize(query_desc.rows); + for (int i = 0; i != idx.rows; ++i) { + auto& row = matches[i]; + row.reserve(k); + for (int j = 0; j != idx.cols; ++j) + row.emplace_back(i, idx.at(i, j), std::sqrt(dst.at(i, j))); + } } diff --git a/src/phg/sfm/homography.cpp b/src/phg/sfm/homography.cpp index 5cbc780..691e9da 100644 --- a/src/phg/sfm/homography.cpp +++ b/src/phg/sfm/homography.cpp @@ -84,8 +84,12 @@ namespace { double w1 = ws1[i]; // 8 elements of matrix + free term as needed by gauss routine -// A.push_back({TODO}); -// A.push_back({TODO}); + A.push_back({.0, .0, .0, + -x0 * w1, -y0 * w1, -w0 * w1, + x0 * y1, y0 * y1, -w0 * y1}); + A.push_back({x0 * w1, y0 * w1, w0 * w1, + .0, .0, .0, + -x0 * x1, -y0 * x1, w0 * x1}); } int res = gauss(A, H); @@ -162,63 +166,62 @@ namespace { throw std::runtime_error("findHomography: points_lhs.size() != points_rhs.size()"); } - // 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(); + 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_samples = 4; + const int n_trials = std::ceil(std::log(1 - 0.999) / std::log(1 - std::pow(0.5, n_samples))); + 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,8 @@ 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 res = T * cv::Vec3d(pt.x, pt.y, 1); + return {res.at(0) / res.at(2), res.at(1) / res.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..c52096e 100644 --- a/src/phg/sfm/panorama_stitcher.cpp +++ b/src/phg/sfm/panorama_stitcher.cpp @@ -4,6 +4,7 @@ #include #include + /* * imgs - список картинок * parent - список индексов, каждый индекс указывает, к какой картинке должна быть приклеена текущая картинка @@ -20,12 +21,30 @@ 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"); + + // здесь надо посчитать вектор Hs + // при этом можно обойтись n_images - 1 вызовами функтора homography_builder + + std::vector used(n_images, false); + const auto dfs = [&](const auto &self, int i) { + if (used[i]) { + return Hs[i]; + } + used[i] = true; + + if (parent[i] == -1) { + Hs[i] = cv::Mat::eye(3, 3, CV_64FC1); + } else { + Hs[i] = homography_builder(imgs[i], imgs[parent[i]]) * self(self, parent[i]); + } + return Hs[i]; + }; + + for (int i = 0; i < n_images; ++i) { + dfs(dfs, i); } + bbox2 bbox; for (int i = 0; i < n_images; ++i) { double w = imgs[i].cols; diff --git a/src/phg/sift/sift.cpp b/src/phg/sift/sift.cpp index 7204771..69b01d0 100755 --- a/src/phg/sift/sift.cpp +++ b/src/phg/sift/sift.cpp @@ -18,50 +18,47 @@ namespace { -cv::Mat upsample2x(const cv::Mat& src) -{ - cv::Mat dst; - cv::resize(src, dst, cv::Size(src.cols * 2, src.rows * 2), 0, 0, cv::INTER_LINEAR); - return dst; -} + cv::Mat upsample2x(const cv::Mat& src) { + cv::Mat dst; + cv::resize(src, dst, cv::Size(src.cols * 2, src.rows * 2), 0, 0, cv::INTER_LINEAR); + return dst; + } -cv::Mat downsample2x(const cv::Mat& src) -{ - int dstW = src.cols / 2; - int dstH = src.rows / 2; - cv::Mat dst(dstH, dstW, src.type()); - const int ch = src.channels(); - - for (int y = 0; y < dstH; y++) { - const float* srcRow = src.ptr(y * 2); - float* dstRow = dst.ptr(y); - for (int x = 0; x < dstW; x++) { - std::copy(srcRow + x * 2 * ch, srcRow + x * 2 * ch + ch, dstRow + x * ch); + cv::Mat downsample2x(const cv::Mat& src) { + int dstW = src.cols / 2; + int dstH = src.rows / 2; + cv::Mat dst(dstH, dstW, src.type()); + const int ch = src.channels(); + + for (int y = 0; y < dstH; y++) { + const float* srcRow = src.ptr(y * 2); + float* dstRow = dst.ptr(y); + for (int x = 0; x < dstW; x++) { + std::copy(srcRow + x * 2 * ch, srcRow + x * 2 * ch + ch, dstRow + x * ch); + } } + + return dst; } - return dst; -} -[[maybe_unused]] cv::Mat downsample2x_avg(const cv::Mat& src) -{ - int dstW = src.cols / 2; - int dstH = src.rows / 2; - cv::Mat dst(dstH, dstW, src.type()); - - for (int y = 0; y < dstH; y++) { - const float* r0 = src.ptr(y * 2); - const float* r1 = src.ptr(y * 2 + 1); - float* dstRow = dst.ptr(y); - for (int x = 0; x < dstW; x++) { - dstRow[x] = (r0[x * 2] + r0[x * 2 + 1] + r1[x * 2] + r1[x * 2 + 1]) * 0.25f; + [[maybe_unused]] cv::Mat downsample2x_avg(const cv::Mat& src) { + int dstW = src.cols / 2; + int dstH = src.rows / 2; + cv::Mat dst(dstH, dstW, src.type()); + + for (int y = 0; y < dstH; y++) { + const float* r0 = src.ptr(y * 2); + const float* r1 = src.ptr(y * 2 + 1); + float* dstRow = dst.ptr(y); + for (int x = 0; x < dstW; x++) { + dstRow[x] = (r0[x * 2] + r0[x * 2 + 1] + r1[x * 2] + r1[x * 2 + 1]) * 0.25f; + } } + return dst; } - return dst; -} } -cv::Mat phg::toGray32F(const cv::Mat& img) -{ +cv::Mat phg::toGray32F(const cv::Mat& img) { cv::Mat gray; if (img.channels() == 3) { cv::cvtColor(img, gray, cv::COLOR_BGR2GRAY); @@ -76,8 +73,7 @@ cv::Mat phg::toGray32F(const cv::Mat& img) return gray_float; } -std::vector phg::buildOctaves(const cv::Mat& img, const phg::SIFTParams& p, int verbose_level) -{ +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; // взятое с потолка значение блюра который уже есть в картинке. используем для того, чтобы не так сильно блюрить базовую картинку и не терять лишний раз фичи @@ -88,14 +84,24 @@ std::vector phg::buildOctaves(const cv::Mat& img, const phg:: // сигнал реального мира: потенциально высочайшего разрешения, можем зумиться почти до молекул, сигма почти нулевая // сигнал с камеры, входное изображение: было произведено усреднение хотя бы по отдельным пикселям матрицы камеры, что соответствует сигме в полпикселя 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 = разность двух) + 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 пикселей, там уже не имеет смысла что-то детектировать + int n_octaves = std::max(1, (int) std::round(std::log2(std::min(img.cols, img.rows))) - + 3); // не даунскейлим дальше размера картинки в 16 пикселей, там уже не имеет смысла что-то детектировать cv::Mat base; - double sigma_base = std::sqrt(sigma0 * sigma0 - sigma_nominal * sigma_nominal); // можно использовать дальше как идею для инкрементального блюра слоев + double sigma_base = std::sqrt(sigma0 * sigma0 - sigma_nominal * + sigma_nominal); // можно использовать дальше как идею для инкрементального блюра слоев cv::GaussianBlur(img, base, cv::Size(), sigma_base, sigma_base); +// cv::Mat img_8u, base_8u; +// img.convertTo(img_8u, CV_8U, 255.0); +// base.convertTo(base_8u, CV_8U, 255.0); +// cv::imwrite("data/debug/test_sift/debug/img.png", img_8u); +// cv::imwrite("data/debug/test_sift/debug/base.png", base_8u); +// cv::imwrite("data/debug/test_sift/debug/diff.png", img_8u-base_8u); + std::vector octaves(n_octaves); for (int o = 0; o < n_octaves; o++) { @@ -108,16 +114,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(даунскейлим текущий слой в два раза, без интерполяции, просто сабсепмлинг); + base = downsample2x(oct.layers[s]); // можно использовать и downsample2x_avg(oct.layers[s]), это позволяет потом заапскейлить слои обратно до оригинального разрешения без сдвига // но потребуется везде изменить формулу для пересчета ключевых точек: pt_upscaled = (pt_downscaled + 0.5) * 2^o - 0.5 @@ -130,29 +136,32 @@ std::vector phg::buildOctaves(const cv::Mat& img, const phg:: return octaves; } -std::vector phg::buildDoG(const std::vector& octaves, const phg::SIFTParams& p, int verbose_level) -{ +std::vector +phg::buildDoG(const std::vector& octaves, const phg::SIFTParams& p, int verbose_level) { std::vector dog(octaves.size()); for (size_t o = 0; o < octaves.size(); o++) { const phg::SIFT::Octave& octave = octaves[o]; dog[o].layers.resize(octave.layers.size() - 1); - // TODO каждый слой дога это разница n+1 и n-й гауссианы + for (size_t i = 0; i < dog[o].layers.size(); i++) { + dog[o].layers[i] = octave.layers[i + 1] - octave.layers[i]; + } } return dog; } -std::vector phg::findScaleSpaceExtrema(const std::vector& dog, const phg::SIFTParams& params, int verbose_level) -{ +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; const double edge_threshold = params.edge_threshold; // чем больше слоев в октаве, тем меньше разница между ними -> компенсируем порог - const float thresh = (float)(contrast_threshold / s); + const float thresh = (float) (contrast_threshold / s); const int border = 5; @@ -162,11 +171,11 @@ std::vector phg::findScaleSpaceExtrema(const std::vector keypoints; - for (int o = 0; o < (int)dog.size(); o++) { + for (int o = 0; o < (int) dog.size(); o++) { int real_octave = o + first_octave; const std::vector& dog_layers = dog[o].layers; - const int n_dog_layers = (int)dog_layers.size(); + const int n_dog_layers = (int) dog_layers.size(); rassert(n_dog_layers == s + 2, 2138971238612312); // итерируемся по внутренним слоям пирамиды, у нас всегда есть предыдущий и следующий сосед @@ -207,17 +216,30 @@ 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,10 +294,10 @@ std::vector phg::findScaleSpaceExtrema(const std::vector phg::findScaleSpaceExtrema(const std::vector ttd_thresh) + break; } // скейлим координаты точек обратно до родных размеров картинки // !!! если выбираем при даунскейле другую политику, с усреднением вместо ресемплинга, то надо здесь применять формулу со сдвигами на полпикселя - float scale = (real_octave >= 0) ? (float)(1 << real_octave) : (1.f / (float)(1 << (-real_octave))); + float scale = (real_octave >= 0) ? (float) (1 << real_octave) : (1.f / (float) (1 + << (-real_octave))); float real_x = (xi + offset[0]) * scale; float real_y = (yi + offset[1]) * scale; float real_layer = li + offset[2]; @@ -303,7 +327,7 @@ std::vector phg::findScaleSpaceExtrema(const std::vector phg::findScaleSpaceExtrema(const std::vector s || xi < border || xi >= cols - border || yi < border || yi >= rows - border) + if (li < 1 || li > s || xi < border || xi >= cols - border || yi < border || + yi >= rows - border) break; } } @@ -340,8 +365,9 @@ std::vector phg::findScaleSpaceExtrema(const std::vector phg::computeOrientations(const std::vector& kpts, const std::vector& octaves, const phg::SIFTParams& params, int verbose_level) -{ +std::vector +phg::computeOrientations(const std::vector& kpts, const std::vector& octaves, + const phg::SIFTParams& params, int verbose_level) { const int s = params.n_octave_layers; const double sigma0 = params.sigma; const int n_bins = params.orient_nbins; @@ -353,24 +379,25 @@ std::vector phg::computeOrientations(const std::vector= 0) ? (float)(1 << real_octave) : (1.f / (float)(1 << (-real_octave))); + float scale = (real_octave >= 0) ? (float) (1 << real_octave) : (1.f / (float) (1 << (-real_octave))); float x = kp.pt.x / scale; float y = kp.pt.y / scale; // найдем радиус ключевой точки в координатах ее октавы - float kp_sigma_octave = (float)(sigma0 * std::pow(2.0, (double)layer / s)); - float sigma_win = 1.5f * kp_sigma_octave; // цитата из lowe: "Each sample added to the histogram is weighted by its gradient magnitude and by a Gaussian-weighted circular window with a σ that is 1.5 times that of the scale of the keypoint." - int radius = (int)std::round(3.f * sigma_win); + float kp_sigma_octave = (float) (sigma0 * std::pow(2.0, (double) layer / s)); + float sigma_win = 1.5f * + kp_sigma_octave; // цитата из lowe: "Each sample added to the histogram is weighted by its gradient magnitude and by a Gaussian-weighted circular window with a σ that is 1.5 times that of the scale of the keypoint." + int radius = (int) std::round(3.f * sigma_win); - int xi = (int)std::round(x); - int yi = (int)std::round(y); + int xi = (int) std::round(x); + int yi = (int) std::round(y); if (xi - radius <= 0 || xi + radius >= img.cols - 1 || yi - radius <= 0 || yi + radius >= img.rows - 1) continue; @@ -379,39 +406,41 @@ 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(gy, gx); + 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 dxf = (float) dx; + float dyf = (float) dy; + float weight = std::exp(-(dxf * dxf + dyf * dyf) / (2.f * sigma_win * sigma_win)); + if (!params.enable_orientation_gaussian_weighting) { + weight = 1.f; + } + + // голосуем в гистограмме направлений. находим два ближайших бина и гладко распределяем голос между ними + // в таком случае, голос попавший близко к границе между бинами, проголосует поровну за оба бина + float bin = angle_deg * n_bins / 360.f; // [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 = 0.f; + } + + histogram[bin0] += mag * weight * (1.f - frac); + histogram[bin1] += mag * weight * frac; } } @@ -436,7 +465,8 @@ std::vector phg::computeOrientations(const std::vector histogram[prev] && histogram[i] > histogram[next] && histogram[i] >= peak_ratio * max_val) { + if (histogram[i] > histogram[prev] && histogram[i] > histogram[next] && + histogram[i] >= peak_ratio * max_val) { float left = histogram[prev]; float center = histogram[i]; float right = histogram[next]; @@ -450,20 +480,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 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 offset = (left - right) / (2.f * (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); } } } @@ -475,8 +505,9 @@ std::vector phg::computeOrientations(const std::vector> phg::computeDescriptors(const std::vector& kpts, const std::vector& octaves, const phg::SIFTParams& params, int verbose_level) -{ +std::pair> +phg::computeDescriptors(const std::vector& kpts, const std::vector& octaves, + const phg::SIFTParams& params, int verbose_level) { const int s = params.n_octave_layers; const double sigma0 = params.sigma; @@ -496,18 +527,18 @@ std::pair> phg::computeDescriptors(const std: const int first_octave = params.upscale_first ? -1 : 0; - for (const cv::KeyPoint& kp : kpts) { + for (const cv::KeyPoint& kp: kpts) { int layer = kp.class_id; int real_octave = kp.octave; int o = real_octave - first_octave; // индекс в массиве octaves const cv::Mat& img = octaves[o].layers[layer]; - float scale = (real_octave >= 0) ? (float)(1 << real_octave) : (1.f / (float)(1 << (-real_octave))); + float scale = (real_octave >= 0) ? (float) (1 << real_octave) : (1.f / (float) (1 << (-real_octave))); float x = kp.pt.x / scale; float y = kp.pt.y / scale; - float kp_sigma_octave = (float)(sigma0 * std::pow(2.0, (double)layer / s)); + float kp_sigma_octave = (float) (sigma0 * std::pow(2.0, (double) layer / s)); // размер патча в котором считаем дескриптор в пикселях октавы float spatial_bin_width = spatial_bin_width_sigmas * kp_sigma_octave; @@ -516,20 +547,20 @@ std::pair> phg::computeDescriptors(const std: // * sqrt(2) для того, чтобы можно было посемплировать патч даже повернутый на 45 градусов ромбиком // * +1 в скобках чтобы можно было семплировать градиенты (а еще зачем?) float half_width = 0.5f * spatial_bin_width * (n_spatial_bins + 1) * std::sqrt(2.f); - int radius = (int)std::round(half_width); + int radius = (int) std::round(half_width); - int xi = (int)std::round(x); - int yi = (int)std::round(y); + int xi = (int) std::round(x); + int yi = (int) std::round(y); if (xi - radius <= 0 || xi + radius >= img.cols - 1 || yi - radius <= 0 || yi + radius >= img.rows - 1) continue; - float kp_angle_rad = kp.angle * (float)CV_PI / 180.f; + float kp_angle_rad = kp.angle * (float) CV_PI / 180.f; float cos_a = std::cos(kp_angle_rad); 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.5f; std::vector desc(n_dims, 0.f); @@ -551,7 +582,7 @@ std::pair> phg::computeDescriptors(const std: float bin_x = rot_x + n_spatial_bins * 0.5f - 0.5f; float bin_y = rot_y + n_spatial_bins * 0.5f - 0.5f; - if (bin_x < -1.f || bin_x >= (float)n_spatial_bins || bin_y < -1.f || bin_y >= (float)n_spatial_bins) + if (bin_x < -1.f || bin_x >= (float) n_spatial_bins || bin_y < -1.f || bin_y >= (float) n_spatial_bins) continue; // градиент (потом все равно будем все нормализовывать, так что можно не нормировать здесь) @@ -564,9 +595,9 @@ std::pair> phg::computeDescriptors(const std: // инвариантность к повороту: повернем направление градиента на угол ключевой точки float angle_invariant = angle - kp_angle_rad; while (angle_invariant < 0.f) - angle_invariant += (float)CV_2PI; - while (angle_invariant >= (float)CV_2PI) - angle_invariant -= (float)CV_2PI; + angle_invariant += (float) CV_2PI; + while (angle_invariant >= (float) CV_2PI) + angle_invariant -= (float) CV_2PI; // подсчет бина в гистограммке градиентов внутри пространственного бина float bin_o = angle_invariant * n_orient_bins / CV_2PI; @@ -574,18 +605,18 @@ 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 по пространственным бинам и по бинам гистограммок трилинейной интерполяцией - int ix0 = (int)std::floor(bin_x); - int iy0 = (int)std::floor(bin_y); - int io0 = (int)std::floor(bin_o); + int ix0 = (int) std::floor(bin_x); + int iy0 = (int) std::floor(bin_y); + int io0 = (int) std::floor(bin_o); float fx = bin_x - ix0; float fy = bin_y - iy0; @@ -609,20 +640,20 @@ 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; } } } } else { - int ix_nearest = (int)std::round(bin_x); - int iy_nearest = (int)std::round(bin_y); - 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 ix_nearest = (int) std::round(bin_x); + int iy_nearest = (int) std::round(bin_y); + 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) { + int idx = (iy_nearest * n_spatial_bins + ix_nearest) * n_orient_bins + io_nearest; + desc[idx] += weighted_mag; } } } @@ -630,22 +661,22 @@ std::pair> phg::computeDescriptors(const std: // нормализуем дескриптор до единичной l2 длины float norm = 0.f; - for (float v : desc) + for (float v: desc) norm += v * v; norm = std::sqrt(norm) + 1e-7f; - for (float& v : desc) + for (float& v: desc) v /= norm; // грохнем слишком большие градиенты и ренормализуем // таким образом один выброс не потянет за собой весь дескриптор и в будущем расстояние с похожим соседом не вырастет сильно - for (float& v : desc) + for (float& v: desc) v = std::min(v, mag_cap); norm = 0.f; - for (float v : desc) + for (float v: desc) norm += v * v; norm = std::sqrt(norm) + 1e-7f; - for (float& v : desc) + for (float& v: desc) v /= norm; if (descriptors.empty()) { @@ -658,14 +689,15 @@ std::pair> phg::computeDescriptors(const std: } if (verbose_level) - std::cout << "descriptors: " << kpts.size() << " -> " << valid_kpts.size() << " keypoints (some discarded due to border)" << std::endl; + std::cout << "descriptors: " << kpts.size() << " -> " << valid_kpts.size() + << " keypoints (some discarded due to border)" << std::endl; - return { descriptors, valid_kpts }; + return {descriptors, valid_kpts}; } -std::vector phg::selectTopKeypoints(const std::vector& kpts, const phg::SIFTParams& params, int verbose_level) -{ - if (params.nfeatures <= 0 || (int)kpts.size() <= params.nfeatures) { +std::vector +phg::selectTopKeypoints(const std::vector& kpts, const phg::SIFTParams& params, int verbose_level) { + if (params.nfeatures <= 0 || (int) kpts.size() <= params.nfeatures) { return kpts; } @@ -673,7 +705,8 @@ std::vector phg::selectTopKeypoints(const std::vector idx(kpts.size()); std::iota(idx.begin(), idx.end(), 0); - std::partial_sort(idx.begin(), idx.begin() + nfeatures, idx.end(), [&kpts](int a, int b) { return std::abs(kpts[a].response) > std::abs(kpts[b].response); }); + std::partial_sort(idx.begin(), idx.begin() + nfeatures, idx.end(), + [&kpts](int a, int b) { return std::abs(kpts[a].response) > std::abs(kpts[b].response); }); idx.resize(nfeatures); std::sort(idx.begin(), idx.end()); @@ -688,9 +721,10 @@ std::vector phg::selectTopKeypoints(const std::vector& kpts, cv::Mat& desc) const -{ - rassert(mask.empty(), 911738571854310); // not implemented, parameter added to match interface with opencv sift implementation +void phg::SIFT::detectAndCompute(const cv::Mat& img, const cv::Mat& mask, std::vector& kpts, + cv::Mat& desc) const { + rassert(mask.empty(), + 911738571854310); // not implemented, parameter added to match interface with opencv sift implementation saveImg("00_input.jpg", img); @@ -701,7 +735,8 @@ void phg::SIFT::detectAndCompute(const cv::Mat& img, const cv::Mat& mask, std::v auto prev_size = gray.size(); gray = upsample2x(gray); if (verbose_level) - std::cout << "upscaled image from " << prev_size.width << "x" << prev_size.height << " to " << gray.cols << "x" << gray.rows << std::endl; + std::cout << "upscaled image from " << prev_size.width << "x" << prev_size.height << " to " << gray.cols + << "x" << gray.rows << std::endl; saveImg("01b_gray_upscaled.png", gray); } @@ -734,8 +769,7 @@ void phg::SIFT::detectAndCompute(const cv::Mat& img, const cv::Mat& mask, std::v // как подкрутить алгоритм, чтобы всегда выдавать ровно запрошенное количество точек (когда это в принципе возможно) но не сильно просесть в производительности? } -void phg::SIFT::saveImg(const std::string& name, const cv::Mat& img) const -{ +void phg::SIFT::saveImg(const std::string& name, const cv::Mat& img) const { if (verbose_level < 2 || debug_folder.empty()) { return; } @@ -749,8 +783,7 @@ void phg::SIFT::saveImg(const std::string& name, const cv::Mat& img) const cv::imwrite(debug_folder + name, out); } -void phg::SIFT::savePyramid(const std::string& name, const std::vector& pyramid, bool normalize) const -{ +void phg::SIFT::savePyramid(const std::string& name, const std::vector& pyramid, bool normalize) const { if (verbose_level < 2 || debug_folder.empty()) { return; } diff --git a/tests/test_matching.cpp b/tests/test_matching.cpp index adaac65..c9e927e 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 diff --git a/tests/test_sift.cpp b/tests/test_sift.cpp index cf3bd7d..f0727b9 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 @@ -278,7 +275,7 @@ void evaluateDetection(const cv::Mat& M, double minRecall, cv::Mat img0 = cv::Ma // что за числа в дескрипторах двух сопоставленных точек, насколько они похожи, // и сверить что расстояние между дескрипторами - это действительно расстояние // между точками в пространстве высокой размерности: -#if 0 +#if 1 if (i % 100 == 0) { #pragma omp critical { @@ -291,6 +288,7 @@ void evaluateDetection(const cv::Mat& M, double minRecall, cv::Mat img0 = cv::Ma std::cout << "sum((d1-d0)^2): " << cv::sum(mul) << std::endl; std::cout << "sqrt(sum((d1-d0)^2)): " << sqrt(cv::sum(mul)[0]) << std::endl; std::cout << "norm: " << cv::norm(d0, d1, cv::NORM_L2) << std::endl; + std::cout << std::endl; } } #endif