From 0c73863f37ace842c3c8251723d683fb23449043 Mon Sep 17 00:00:00 2001 From: SheLesTT Date: Thu, 5 Mar 2026 17:22:31 +0300 Subject: [PATCH 1/5] tmp --- src/phg/sift/sift.cpp | 186 +++++++++++------- tests/test_sift.cpp | 446 +++++++++++++++++++++--------------------- 2 files changed, 336 insertions(+), 296 deletions(-) diff --git a/src/phg/sift/sift.cpp b/src/phg/sift/sift.cpp index 72047711..2606641c 100755 --- a/src/phg/sift/sift.cpp +++ b/src/phg/sift/sift.cpp @@ -86,7 +86,7 @@ std::vector phg::buildOctaves(const cv::Mat& img, const phg:: // when sampling a continuous signal into a discrete grid. It is a mathematically grounded assumption for digital images. // общая идея в том, что у нас есть какой-то сигнал реального мира, и есть входное изображение // сигнал реального мира: потенциально высочайшего разрешения, можем зумиться почти до молекул, сигма почти нулевая - // сигнал с камеры, входное изображение: было произведено усреднение хотя бы по отдельным пикселям матрицы камеры, что соответствует сигме в полпикселя + // сигнал с камеры, входное изображение: было произведено усреднение хотя бы по отдельным пикселям матрицы камеры, что соответствует сигме в полпикс/еля 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 = разность двух) @@ -107,17 +107,19 @@ std::vector phg::buildOctaves(const cv::Mat& img, const phg:: // для простоты в каждой октаве будем каждый раз блюрить базовую картинку с полной сигмой // можно подумать, как сделать эффективнее - для построения n+1 слоя доблюревать уже поблюренный n-ый слой, так чтобы в итоге получилась такая же сигма // это будет немного быстрее, тк нужно более маленькое ядро свертки на каждый шаг + double sigma_cur = sigma0; for (int i = 1; i < n_layers; i++) { - // TODO double sigma_layer = sigma0 * корень из двух нужной степени, чтобы при i==s получали удвоение базового блюра; - // // вычтем sigma0 чтобы размыть ровно до нужной суммарной сигмы + double sigma_layer_full = sigma_cur * pow(2.0, 1.0/s); + double sigma_layer = std::sqrt(sigma_layer_full * sigma_layer_full - sigma_cur * sigma_cur); // можно использовать дальше как идею для инкрементального блюра слоев // TODO sigma_layer = ... (вычитаем как в sigma base); - // cv::GaussianBlur(oct.layers[0], oct.layers[i], cv::Size(), sigma_layer, sigma_layer); + cv::GaussianBlur(oct.layers[0], oct.layers[i], cv::Size(), sigma_layer, sigma_layer); + sigma_cur = sigma_layer; } // подготавливаем базовый слой для следующей октавы 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 @@ -138,7 +140,9 @@ std::vector phg::buildDoG(const std::vector phg::findScaleSpaceExtrema(const std::vector phg::findScaleSpaceExtrema(const std::vector(yi, xi) - pL.at(yi, xi)) * 0.5f; // гессиан - float dxx, dxy, dyy, dxs, dys, dss; -// float dxx = cL.at(yi, xi + 1) + cL.at(yi, xi - 1) - 2.f * resp_center; -// float dyy = TODO; -// float dss = TODO; -// -// float dxy = (cL.at(yi + 1, xi + 1) - cL.at(yi + 1, xi - 1) - cL.at(yi - 1, xi + 1) + cL.at(yi - 1, xi - 1)) * 0.25f; -// float dxs = TODO; -// float dys = TODO; + float dxx = cL.at(yi, xi + 1) + cL.at(yi, xi - 1) - 2.f * resp_center; + float dyy = cL.at(yi + 1, xi) + cL.at(yi - 1, xi) - 2.f * resp_center; + float dss = nL.at(yi, xi) + pL.at(yi, xi) - 2.f * resp_center; + + float dxy = (cL.at(yi + 1, xi + 1) - cL.at(yi + 1, xi - 1) - cL.at(yi - 1, xi + 1) + cL.at(yi - 1, xi - 1)) * 0.25f; + float dxs = (nL.at(yi, xi + 1) - nL.at(yi, xi - 1) - pL.at(yi, xi + 1) + pL.at(yi, xi - 1)) * 0.25f; + float dys = (nL.at(yi + 1, xi) - nL.at(yi - 1, xi) - pL.at(yi + 1, xi) + pL.at(yi - 1, xi)) * 0.25f; cv::Matx33f H(dxx, dxy, dxs, dxy, dyy, dys, dxs, dys, dss); @@ -273,8 +305,8 @@ std::vector phg::findScaleSpaceExtrema(const std::vector phg::findScaleSpaceExtrema(const std::vector pow((r+1),2)/r){ + break; + } } // скейлим координаты точек обратно до родных размеров картинки @@ -379,39 +412,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; -// } + 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] // -// histogram[bin0] += TODO; -// histogram[bin1] += TODO; + 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; + if (bin >= n_bins) bin -= n_bins; + int bin0 = (int) bin; + int bin1 = (bin0 + 1) % n_bins; + + float frac = bin - bin0; + if (!params.enable_orientation_bin_interpolation) { + frac = 0.f; + } + + histogram[bin0] += (1.0f - frac) * mag * weight; + histogram[bin1] += frac * mag * weight; } } @@ -450,20 +483,22 @@ std::vector phg::computeOrientations(const std::vector a = (left + right - 2 * center) / 2 // f(1) - f(-1) = 2b -> b = (right - left) / 2 -// float offset = TODO; -// if (!params.enable_orientation_subpixel_localization) { -// offset = 0.f; -// } -// -// float bin_real = i + offset; -// if (bin_real < 0.f) bin_real += n_bins; -// if (bin_real >= n_bins) bin_real -= n_bins; -// -// float angle = bin_real * 360.f / n_bins; -// -// cv::KeyPoint new_kp = kp; -// new_kp.angle = angle; -// oriented_kpts.push_back(new_kp); + float a = (left + right - 2.f * center) / 2.f; + float b = (right - left) / 2.f; + float offset = -b / (2.f * a); + if (!params.enable_orientation_subpixel_localization) { + offset = 0.f; + } + + float bin_real = i + offset; + if (bin_real < 0.f) bin_real += n_bins; + if (bin_real >= n_bins) bin_real -= n_bins; + + float angle = bin_real * 360.f / n_bins; + + cv::KeyPoint new_kp = kp; + new_kp.angle = angle; + oriented_kpts.push_back(new_kp); } } } @@ -574,11 +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;; } } } @@ -620,9 +655,8 @@ std::pair> phg::computeDescriptors(const std: int io_nearest = (int)std::round(bin_o) % n_orient_bins; if (ix_nearest >= 0 && ix_nearest < n_spatial_bins && iy_nearest >= 0 && iy_nearest < n_spatial_bins) { - // TODO uncomment -// int idx = (iy_nearest * n_spatial_bins + ix_nearest) * n_orient_bins + io_nearest; -// desc[idx] += weighted_mag; + int idx = (iy_nearest * n_spatial_bins + ix_nearest) * n_orient_bins + io_nearest; + desc[dx] += weighted_mag; } } } diff --git a/tests/test_sift.cpp b/tests/test_sift.cpp index cf3bd7df..ed393a95 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 @@ -139,8 +139,10 @@ double diffAngles(double angle0, double angle1) // На вход передается матрица описывающая преобразование картинки (сдвиг, поворот, масштабирование или их комбинация), допустимый процент Recall, и опционально можно тестировать другую картинку void evaluateDetection(const cv::Mat& M, double minRecall, cv::Mat img0 = cv::Mat()) { + #include + std::cout << std::filesystem::current_path() << std::endl; if (img0.empty()) { - img0 = cv::imread("data/src/test_sift/unicorn.png"); // грузим картинку по умолчанию + img0 = cv::imread("/home/denisden/photo/PhotogrammetryTasks2026/data/src/test_sift/unicorn.png"); // грузим картинку по умолчанию } ASSERT_FALSE(img0.empty()); // проверка что картинка была загружена @@ -212,7 +214,11 @@ void evaluateDetection(const cv::Mat& M, double minRecall, cv::Mat img0 = cv::Ma for (size_t i = 0; i < kps0.size(); ++i) { ps0[i] = kps0[i].pt; } + + std::cout << "before" <Working directory) указывает на корневую папку проекта - - double angleDegreesClockwise = -40; - double scale = 1.0; - double minRecall = 0.75; - evaluateDetection(cv::getRotationMatrix2D(cv::Point(jesu19.cols / 2, jesu19.rows / 2), -angleDegreesClockwise, scale), minRecall, jesu19); -} - -TEST(SIFT, DetectionSmokeTest) -{ -#if ENABLE_MY_SIFT_TESTING - phg::SIFTParams p; - phg::SIFT sift(p, 2, "data/debug/test_sift/debug/"); - - cv::Mat img = cv::imread("data/src/test_sift/mysh1.jpg"); - cv::resize(img, img, img.size() / 4, 0, 0, cv::INTER_AREA); - - std::vector kpts; - cv::Mat desc; - sift.detectAndCompute(img, kpts, desc); -#else - std::cout << "ENABLE_MY_SIFT_TESTING is disabled, test skipped" << std::endl; -#endif -} +// TEST(SIFT, MovedTheSameImage) +// { +// double minRecall = 0.75; +// evaluateDetection(createTranslationMatrix(0.0, 0.0), minRecall); +// } + +// TEST(SIFT, MovedImageRight) +// { +// double minRecall = 0.75; +// evaluateDetection(createTranslationMatrix(50.0, 0.0), minRecall); +// } + +// TEST(SIFT, MovedImageLeft) +// { +// double minRecall = 0.75; +// evaluateDetection(createTranslationMatrix(-50.0, 0.0), minRecall); +// } + +// TEST(SIFT, MovedImageUpHalfPixel) +// { +// double minRecall = 0.75; +// evaluateDetection(createTranslationMatrix(0.0, -50.5), minRecall); +// } + +// TEST(SIFT, MovedImageDownHalfPixel) +// { +// double minRecall = 0.75; +// evaluateDetection(createTranslationMatrix(0.0, 50.5), minRecall); +// } + +// TEST(SIFT, Rotate10) +// { +// double angleDegreesClockwise = 10; +// double scale = 1.0; +// double minRecall = 0.60; +// evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +// } + +// TEST(SIFT, Rotate20) +// { +// double angleDegreesClockwise = 20; +// double scale = 1.0; +// double minRecall = 0.60; +// evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +// } + +// TEST(SIFT, Rotate30) +// { +// double angleDegreesClockwise = 30; +// double scale = 1.0; +// double minRecall = 0.60; +// evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +// } + +// TEST(SIFT, Rotate40) +// { +// double angleDegreesClockwise = 40; +// double scale = 1.0; +// double minRecall = 0.60; +// evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +// } + +// TEST(SIFT, Rotate45) +// { +// double angleDegreesClockwise = 45; +// double scale = 1.0; +// double minRecall = 0.60; +// evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +// } + +// TEST(SIFT, Rotate90) +// { +// double angleDegreesClockwise = 90; +// double scale = 1.0; +// double minRecall = 0.75; +// evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +// } + +// TEST(SIFT, Scale50) +// { +// double angleDegreesClockwise = 0; +// double scale = 0.5; +// double minRecall = 0.40; +// evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +// } + +// TEST(SIFT, Scale70) +// { +// double angleDegreesClockwise = 0; +// double scale = 0.7; +// double minRecall = 0.40; +// evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +// } + +// TEST(SIFT, Scale90) +// { +// double angleDegreesClockwise = 0; +// double scale = 0.9; +// double minRecall = 0.60; +// evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +// } + +// TEST(SIFT, Scale110) +// { +// double angleDegreesClockwise = 0; +// double scale = 1.1; +// double minRecall = 0.60; +// evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +// } + +// TEST(SIFT, Scale130) +// { +// double angleDegreesClockwise = 0; +// double scale = 1.3; +// double minRecall = 0.50; +// evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +// } + +// TEST(SIFT, Scale150) +// { +// double angleDegreesClockwise = 0; +// double scale = 1.5; +// double minRecall = 0.50; +// evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +// } + +// TEST(SIFT, Scale175) +// { +// double angleDegreesClockwise = 0; +// double scale = 1.75; +// double minRecall = 0.40; +// evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +// } + +// TEST(SIFT, Scale200) +// { +// double angleDegreesClockwise = 0; +// double scale = 2.0; +// double minRecall = 0.20; +// evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +// } + +// TEST(SIFT, Rotate10Scale90) +// { +// double angleDegreesClockwise = 10; +// double scale = 0.9; +// double minRecall = 0.65; +// evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +// } + +// TEST(SIFT, Rotate30Scale75) +// { +// double angleDegreesClockwise = 30; +// double scale = 0.75; +// double minRecall = 0.50; +// evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +// } + +// TEST(SIFT, HerzJesu19RotateM40) +// { +// cv::Mat jesu19 = cv::imread("data/src/test_sift/herzjesu19.png"); + +// ASSERT_FALSE(jesu19.empty()); // проверка что картинка была загружена +// // убедитесь что рабочая папка (Edit Configurations...->Working directory) указывает на корневую папку проекта + +// double angleDegreesClockwise = -40; +// double scale = 1.0; +// double minRecall = 0.75; +// evaluateDetection(cv::getRotationMatrix2D(cv::Point(jesu19.cols / 2, jesu19.rows / 2), -angleDegreesClockwise, scale), minRecall, jesu19); +// } + +// TEST(SIFT, DetectionSmokeTest) +// { +// #if ENABLE_MY_SIFT_TESTING +// phg::SIFTParams p; +// phg::SIFT sift(p, 2, "data/debug/test_sift/debug/"); + +// cv::Mat img = cv::imread("data/src/test_sift/mysh1.jpg"); +// cv::resize(img, img, img.size() / 4, 0, 0, cv::INTER_AREA); + +// std::vector kpts; +// cv::Mat desc; +// sift.detectAndCompute(img, kpts, desc); +// #else +// std::cout << "ENABLE_MY_SIFT_TESTING is disabled, test skipped" << std::endl; +// #endif +// } namespace fs = std::filesystem; @@ -899,42 +905,42 @@ TEST(SIFT, DetectionDescriptionSteps) #endif } -TEST(SIFT, PairMatching) -{ -#if ENABLE_MY_SIFT_TESTING - cv::Mat img1 = cv::imread("data/src/test_sift/mysh2.jpg"); - ASSERT_FALSE(img1.empty()); +// TEST(SIFT, PairMatching) +// { +// #if ENABLE_MY_SIFT_TESTING +// cv::Mat img1 = cv::imread("data/src/test_sift/mysh2.jpg"); +// ASSERT_FALSE(img1.empty()); - cv::Mat img2 = cv::imread("data/src/test_sift/mysh3.jpg"); - ASSERT_FALSE(img2.empty()); +// cv::Mat img2 = cv::imread("data/src/test_sift/mysh3.jpg"); +// ASSERT_FALSE(img2.empty()); - cv::resize(img1, img1, img1.size() / 2, 0, 0, cv::INTER_AREA); - cv::resize(img2, img2, img2.size() / 2, 0, 0, cv::INTER_AREA); - std::cout << "image sizes: " << img1.size() << ", " << img2.size() << std::endl; +// cv::resize(img1, img1, img1.size() / 2, 0, 0, cv::INTER_AREA); +// cv::resize(img2, img2, img2.size() / 2, 0, 0, cv::INTER_AREA); +// std::cout << "image sizes: " << img1.size() << ", " << img2.size() << std::endl; - phg::SIFTParams params; - params.nfeatures = 10000; +// phg::SIFTParams params; +// params.nfeatures = 10000; - std::cout << "matching using opencv orb..." << std::endl; - auto orb_cv = cv::ORB::create(params.nfeatures); - evaluateMatching(*orb_cv, img1, img2, "data/debug/test_sift/SIFT/Matches_ORB.jpg"); +// std::cout << "matching using opencv orb..." << std::endl; +// auto orb_cv = cv::ORB::create(params.nfeatures); +// evaluateMatching(*orb_cv, img1, img2, "data/debug/test_sift/SIFT/Matches_ORB.jpg"); - std::cout << "matching using opencv sift..." << std::endl; - auto sift_cv = cv::SIFT::create(params.nfeatures, params.n_octave_layers, params.contrast_threshold, params.edge_threshold); - MatchingPairData data_cv = evaluateMatching(*sift_cv, img1, img2, "data/debug/test_sift/SIFT/Matches_SIFTOCV.jpg"); +// std::cout << "matching using opencv sift..." << std::endl; +// auto sift_cv = cv::SIFT::create(params.nfeatures, params.n_octave_layers, params.contrast_threshold, params.edge_threshold); +// MatchingPairData data_cv = evaluateMatching(*sift_cv, img1, img2, "data/debug/test_sift/SIFT/Matches_SIFTOCV.jpg"); - std::cout << "matching using my sift..." << std::endl; - phg::SIFT sift(params); - MatchingPairData data = evaluateMatching(sift, img1, img2, "data/debug/test_sift/SIFT/Matches_SIFT_MY.jpg"); +// std::cout << "matching using my sift..." << std::endl; +// phg::SIFT sift(params); +// MatchingPairData data = evaluateMatching(sift, img1, img2, "data/debug/test_sift/SIFT/Matches_SIFT_MY.jpg"); - double thresh = 0.8; // expect at least 80% of opencv sift points & matches - EXPECT_GE(data.npoints1, thresh * data_cv.npoints1); - EXPECT_GE(data.npoints2, thresh * data_cv.npoints2); - EXPECT_GE(data.nmatches, thresh * data_cv.nmatches); +// double thresh = 0.8; // expect at least 80% of opencv sift points & matches +// EXPECT_GE(data.npoints1, thresh * data_cv.npoints1); +// EXPECT_GE(data.npoints2, thresh * data_cv.npoints2); +// EXPECT_GE(data.nmatches, thresh * data_cv.nmatches); - std::cout << "Final score: " << data.nmatches << std::endl; -#else - std::cout << "ENABLE_MY_SIFT_TESTING is disabled, test skipped" << std::endl; - std::cout << "Final score: UNKNOWN" << std::endl; -#endif -} +// std::cout << "Final score: " << data.nmatches << std::endl; +// #else +// std::cout << "ENABLE_MY_SIFT_TESTING is disabled, test skipped" << std::endl; +// std::cout << "Final score: UNKNOWN" << std::endl; +// #endif +// } From 6ebdf0d505394826747574ebd6d1df78301777b8 Mon Sep 17 00:00:00 2001 From: SheLesTT Date: Fri, 6 Mar 2026 08:26:13 +0300 Subject: [PATCH 2/5] sift --- src/phg/sift/sift.cpp | 16 +- tests/test_sift.cpp | 439 +++++++++++++++++++++--------------------- 2 files changed, 227 insertions(+), 228 deletions(-) diff --git a/src/phg/sift/sift.cpp b/src/phg/sift/sift.cpp index 2606641c..5e90fef3 100755 --- a/src/phg/sift/sift.cpp +++ b/src/phg/sift/sift.cpp @@ -109,11 +109,13 @@ std::vector phg::buildOctaves(const cv::Mat& img, const phg:: // это будет немного быстрее, тк нужно более маленькое ядро свертки на каждый шаг double sigma_cur = sigma0; for (int i = 1; i < n_layers; i++) { - double sigma_layer_full = sigma_cur * pow(2.0, 1.0/s); + double k = pow(2.0, 1.0/s); + double sigma_layer_full = sigma_cur *k ; double sigma_layer = std::sqrt(sigma_layer_full * sigma_layer_full - sigma_cur * sigma_cur); // можно использовать дальше как идею для инкрементального блюра слоев // TODO sigma_layer = ... (вычитаем как в sigma base); - cv::GaussianBlur(oct.layers[0], oct.layers[i], cv::Size(), sigma_layer, sigma_layer); - sigma_cur = sigma_layer; + oct.layers[i] = base.clone(); + cv::GaussianBlur(oct.layers[i-1], oct.layers[i], cv::Size(), sigma_layer, sigma_layer); + sigma_cur = sigma_layer_full; } // подготавливаем базовый слой для следующей октавы @@ -141,7 +143,7 @@ std::vector phg::buildDoG(const std::vector phg::findScaleSpaceExtrema(const std::vector phg::computeOrientations(const std::vectorWorking directory) указывает на корневую папку проекта - -// double angleDegreesClockwise = -40; -// double scale = 1.0; -// double minRecall = 0.75; -// evaluateDetection(cv::getRotationMatrix2D(cv::Point(jesu19.cols / 2, jesu19.rows / 2), -angleDegreesClockwise, scale), minRecall, jesu19); -// } - -// TEST(SIFT, DetectionSmokeTest) -// { -// #if ENABLE_MY_SIFT_TESTING -// phg::SIFTParams p; -// phg::SIFT sift(p, 2, "data/debug/test_sift/debug/"); - -// cv::Mat img = cv::imread("data/src/test_sift/mysh1.jpg"); -// cv::resize(img, img, img.size() / 4, 0, 0, cv::INTER_AREA); - -// std::vector kpts; -// cv::Mat desc; -// sift.detectAndCompute(img, kpts, desc); -// #else -// std::cout << "ENABLE_MY_SIFT_TESTING is disabled, test skipped" << std::endl; -// #endif -// } +TEST(SIFT, MovedTheSameImage) +{ + double minRecall = 0.75; + evaluateDetection(createTranslationMatrix(0.0, 0.0), minRecall); +} + +TEST(SIFT, MovedImageRight) +{ + double minRecall = 0.75; + evaluateDetection(createTranslationMatrix(50.0, 0.0), minRecall); +} + +TEST(SIFT, MovedImageLeft) +{ + double minRecall = 0.75; + evaluateDetection(createTranslationMatrix(-50.0, 0.0), minRecall); +} + +TEST(SIFT, MovedImageUpHalfPixel) +{ + double minRecall = 0.75; + evaluateDetection(createTranslationMatrix(0.0, -50.5), minRecall); +} + +TEST(SIFT, MovedImageDownHalfPixel) +{ + double minRecall = 0.75; + evaluateDetection(createTranslationMatrix(0.0, 50.5), minRecall); +} + +TEST(SIFT, Rotate10) +{ + double angleDegreesClockwise = 10; + double scale = 1.0; + double minRecall = 0.60; + evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +} + +TEST(SIFT, Rotate20) +{ + double angleDegreesClockwise = 20; + double scale = 1.0; + double minRecall = 0.60; + evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +} + +TEST(SIFT, Rotate30) +{ + double angleDegreesClockwise = 30; + double scale = 1.0; + double minRecall = 0.60; + evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +} + +TEST(SIFT, Rotate40) +{ + double angleDegreesClockwise = 40; + double scale = 1.0; + double minRecall = 0.60; + evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +} + +TEST(SIFT, Rotate45) +{ + double angleDegreesClockwise = 45; + double scale = 1.0; + double minRecall = 0.60; + evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +} + +TEST(SIFT, Rotate90) +{ + double angleDegreesClockwise = 90; + double scale = 1.0; + double minRecall = 0.75; + evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +} + +TEST(SIFT, Scale50) +{ + double angleDegreesClockwise = 0; + double scale = 0.5; + double minRecall = 0.40; + evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +} + +TEST(SIFT, Scale70) +{ + double angleDegreesClockwise = 0; + double scale = 0.7; + double minRecall = 0.40; + evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +} + +TEST(SIFT, Scale90) +{ + double angleDegreesClockwise = 0; + double scale = 0.9; + double minRecall = 0.60; + evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +} + +TEST(SIFT, Scale110) +{ + double angleDegreesClockwise = 0; + double scale = 1.1; + double minRecall = 0.60; + evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +} + +TEST(SIFT, Scale130) +{ + double angleDegreesClockwise = 0; + double scale = 1.3; + double minRecall = 0.50; + evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +} + +TEST(SIFT, Scale150) +{ + double angleDegreesClockwise = 0; + double scale = 1.5; + double minRecall = 0.50; + evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +} + +TEST(SIFT, Scale175) +{ + double angleDegreesClockwise = 0; + double scale = 1.75; + double minRecall = 0.40; + evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +} + +TEST(SIFT, Scale200) +{ + double angleDegreesClockwise = 0; + double scale = 2.0; + double minRecall = 0.20; + evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +} + +TEST(SIFT, Rotate10Scale90) +{ + double angleDegreesClockwise = 10; + double scale = 0.9; + double minRecall = 0.65; + evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +} + +TEST(SIFT, Rotate30Scale75) +{ + double angleDegreesClockwise = 30; + double scale = 0.75; + double minRecall = 0.50; + evaluateDetection(cv::getRotationMatrix2D(cv::Point(200, 256), -angleDegreesClockwise, scale), minRecall); +} + +TEST(SIFT, HerzJesu19RotateM40) +{ + cv::Mat jesu19 = cv::imread("data/src/test_sift/herzjesu19.png"); + + ASSERT_FALSE(jesu19.empty()); // проверка что картинка была загружена + // убедитесь что рабочая папка (Edit Configurations...->Working directory) указывает на корневую папку проекта + + double angleDegreesClockwise = -40; + double scale = 1.0; + double minRecall = 0.75; + evaluateDetection(cv::getRotationMatrix2D(cv::Point(jesu19.cols / 2, jesu19.rows / 2), -angleDegreesClockwise, scale), minRecall, jesu19); +} + +TEST(SIFT, DetectionSmokeTest) +{ +#if ENABLE_MY_SIFT_TESTING + phg::SIFTParams p; + phg::SIFT sift(p, 2, "data/debug/test_sift/debug/"); + + cv::Mat img = cv::imread("data/src/test_sift/mysh1.jpg"); + cv::resize(img, img, img.size() / 4, 0, 0, cv::INTER_AREA); + + std::vector kpts; + cv::Mat desc; + sift.detectAndCompute(img, kpts, desc); +#else + std::cout << "ENABLE_MY_SIFT_TESTING is disabled, test skipped" << std::endl; +#endif +} namespace fs = std::filesystem; @@ -905,42 +902,42 @@ TEST(SIFT, DetectionDescriptionSteps) #endif } -// TEST(SIFT, PairMatching) -// { -// #if ENABLE_MY_SIFT_TESTING -// cv::Mat img1 = cv::imread("data/src/test_sift/mysh2.jpg"); -// ASSERT_FALSE(img1.empty()); +TEST(SIFT, PairMatching) +{ +#if ENABLE_MY_SIFT_TESTING + cv::Mat img1 = cv::imread("data/src/test_sift/mysh2.jpg"); + ASSERT_FALSE(img1.empty()); -// cv::Mat img2 = cv::imread("data/src/test_sift/mysh3.jpg"); -// ASSERT_FALSE(img2.empty()); + cv::Mat img2 = cv::imread("data/src/test_sift/mysh3.jpg"); + ASSERT_FALSE(img2.empty()); -// cv::resize(img1, img1, img1.size() / 2, 0, 0, cv::INTER_AREA); -// cv::resize(img2, img2, img2.size() / 2, 0, 0, cv::INTER_AREA); -// std::cout << "image sizes: " << img1.size() << ", " << img2.size() << std::endl; + cv::resize(img1, img1, img1.size() / 2, 0, 0, cv::INTER_AREA); + cv::resize(img2, img2, img2.size() / 2, 0, 0, cv::INTER_AREA); + std::cout << "image sizes: " << img1.size() << ", " << img2.size() << std::endl; -// phg::SIFTParams params; -// params.nfeatures = 10000; + phg::SIFTParams params; + params.nfeatures = 10000; -// std::cout << "matching using opencv orb..." << std::endl; -// auto orb_cv = cv::ORB::create(params.nfeatures); -// evaluateMatching(*orb_cv, img1, img2, "data/debug/test_sift/SIFT/Matches_ORB.jpg"); + std::cout << "matching using opencv orb..." << std::endl; + auto orb_cv = cv::ORB::create(params.nfeatures); + evaluateMatching(*orb_cv, img1, img2, "data/debug/test_sift/SIFT/Matches_ORB.jpg"); -// std::cout << "matching using opencv sift..." << std::endl; -// auto sift_cv = cv::SIFT::create(params.nfeatures, params.n_octave_layers, params.contrast_threshold, params.edge_threshold); -// MatchingPairData data_cv = evaluateMatching(*sift_cv, img1, img2, "data/debug/test_sift/SIFT/Matches_SIFTOCV.jpg"); + std::cout << "matching using opencv sift..." << std::endl; + auto sift_cv = cv::SIFT::create(params.nfeatures, params.n_octave_layers, params.contrast_threshold, params.edge_threshold); + MatchingPairData data_cv = evaluateMatching(*sift_cv, img1, img2, "data/debug/test_sift/SIFT/Matches_SIFTOCV.jpg"); -// std::cout << "matching using my sift..." << std::endl; -// phg::SIFT sift(params); -// MatchingPairData data = evaluateMatching(sift, img1, img2, "data/debug/test_sift/SIFT/Matches_SIFT_MY.jpg"); + std::cout << "matching using my sift..." << std::endl; + phg::SIFT sift(params); + MatchingPairData data = evaluateMatching(sift, img1, img2, "data/debug/test_sift/SIFT/Matches_SIFT_MY.jpg"); -// double thresh = 0.8; // expect at least 80% of opencv sift points & matches -// EXPECT_GE(data.npoints1, thresh * data_cv.npoints1); -// EXPECT_GE(data.npoints2, thresh * data_cv.npoints2); -// EXPECT_GE(data.nmatches, thresh * data_cv.nmatches); + double thresh = 0.8; // expect at least 80% of opencv sift points & matches + EXPECT_GE(data.npoints1, thresh * data_cv.npoints1); + EXPECT_GE(data.npoints2, thresh * data_cv.npoints2); + EXPECT_GE(data.nmatches, thresh * data_cv.nmatches); -// std::cout << "Final score: " << data.nmatches << std::endl; -// #else -// std::cout << "ENABLE_MY_SIFT_TESTING is disabled, test skipped" << std::endl; -// std::cout << "Final score: UNKNOWN" << std::endl; -// #endif -// } + std::cout << "Final score: " << data.nmatches << std::endl; +#else + std::cout << "ENABLE_MY_SIFT_TESTING is disabled, test skipped" << std::endl; + std::cout << "Final score: UNKNOWN" << std::endl; +#endif +} From a2362ce524b5f40896b5b1fc1dd311de092a2e5c Mon Sep 17 00:00:00 2001 From: SheLesTT Date: Fri, 6 Mar 2026 08:48:23 +0300 Subject: [PATCH 3/5] fix path --- tests/test_sift.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_sift.cpp b/tests/test_sift.cpp index 9a8a5c48..290ed831 100755 --- a/tests/test_sift.cpp +++ b/tests/test_sift.cpp @@ -142,7 +142,7 @@ void evaluateDetection(const cv::Mat& M, double minRecall, cv::Mat img0 = cv::Ma #include std::cout << std::filesystem::current_path() << std::endl; if (img0.empty()) { - img0 = cv::imread("/home/denisden/photo/PhotogrammetryTasks2026/data/src/test_sift/unicorn.png"); // грузим картинку по умолчанию + img0 = cv::imread("/data/src/test_sift/unicorn.png"); // грузим картинку по умолчанию } ASSERT_FALSE(img0.empty()); // проверка что картинка была загружена From 4c3c75b9c6c801e8c0833d342036d24a2cbfb8a7 Mon Sep 17 00:00:00 2001 From: SheLesTT Date: Fri, 6 Mar 2026 08:59:09 +0300 Subject: [PATCH 4/5] fix path --- tests/test_sift.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_sift.cpp b/tests/test_sift.cpp index 290ed831..b347faa0 100755 --- a/tests/test_sift.cpp +++ b/tests/test_sift.cpp @@ -142,7 +142,7 @@ void evaluateDetection(const cv::Mat& M, double minRecall, cv::Mat img0 = cv::Ma #include std::cout << std::filesystem::current_path() << std::endl; if (img0.empty()) { - img0 = cv::imread("/data/src/test_sift/unicorn.png"); // грузим картинку по умолчанию + img0 = cv::imread("data/src/test_sift/unicorn.png"); // грузим картинку по умолчанию } ASSERT_FALSE(img0.empty()); // проверка что картинка была загружена From 318361efa7920337d59c78b23ad270d0ad11de0e Mon Sep 17 00:00:00 2001 From: SheLesTT Date: Wed, 18 Mar 2026 16:07:58 +0300 Subject: [PATCH 5/5] task02 --- src/phg/matching/descriptor_matcher.cpp | 107 ++++++++++++-------- src/phg/matching/flann_matcher.cpp | 17 +++- src/phg/sfm/homography.cpp | 125 ++++++++++++++---------- src/phg/sfm/panorama_stitcher.cpp | 31 +++++- 4 files changed, 178 insertions(+), 102 deletions(-) diff --git a/src/phg/matching/descriptor_matcher.cpp b/src/phg/matching/descriptor_matcher.cpp index f4bcd871..6db87bca 100644 --- a/src/phg/matching/descriptor_matcher.cpp +++ b/src/phg/matching/descriptor_matcher.cpp @@ -1,14 +1,20 @@ #include "descriptor_matcher.h" - +#include #include #include "flann_factory.h" - +using namespace std; +#include 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(auto &match : matches) { + if(match[0].distance < 0.7 * match[1].distance){ + filtered_matches.push_back(match[0]); + } + } } @@ -35,42 +41,61 @@ void phg::DescriptorMatcher::filterMatchesClusters(const std::vector points_query.at(i) = keypoints_query[matches[i].queryIdx].pt; points_train.at(i) = keypoints_train[matches[i].trainIdx].pt; } -// -// // размерность всего 2, так что точное KD-дерево -// std::shared_ptr index_params = flannKdTreeIndexParams(TODO); -// std::shared_ptr search_params = flannKsTreeSearchParams(TODO); -// -// std::shared_ptr index_query = flannKdTreeIndex(points_query, index_params); -// std::shared_ptr index_train = flannKdTreeIndex(points_train, index_params); -// -// // для каждой точки найти total neighbors ближайших соседей -// cv::Mat indices_query(n_matches, total_neighbours, CV_32SC1); -// cv::Mat distances2_query(n_matches, total_neighbours, CV_32FC1); -// cv::Mat indices_train(n_matches, total_neighbours, CV_32SC1); -// cv::Mat distances2_train(n_matches, total_neighbours, CV_32FC1); -// -// index_query->knnSearch(points_query, indices_query, distances2_query, total_neighbours, *search_params); -// index_train->knnSearch(points_train, indices_train, distances2_train, total_neighbours, *search_params); -// -// // оценить радиус поиска для каждой картинки -// // NB: radius2_query, radius2_train: квадраты радиуса! -// float radius2_query, radius2_train; -// { -// std::vector max_dists2_query(n_matches); -// std::vector max_dists2_train(n_matches); -// for (int i = 0; i < n_matches; ++i) { -// max_dists2_query[i] = distances2_query.at(i, total_neighbours - 1); -// max_dists2_train[i] = distances2_train.at(i, total_neighbours - 1); -// } -// -// int median_pos = n_matches / 2; -// std::nth_element(max_dists2_query.begin(), max_dists2_query.begin() + median_pos, max_dists2_query.end()); -// std::nth_element(max_dists2_train.begin(), max_dists2_train.begin() + median_pos, max_dists2_train.end()); -// -// radius2_query = max_dists2_query[median_pos] * radius_limit_scale * radius_limit_scale; -// radius2_train = max_dists2_train[median_pos] * radius_limit_scale * radius_limit_scale; -// } -// -// метч остается, если левое и правое множества первых total_neighbors соседей в радиусах поиска(radius2_query, radius2_train) имеют как минимум consistent_matches общих элементов -// // TODO заполнить filtered_matches + + // размерность всего 2, так что точное KD-дерево + std::shared_ptr index_params = flannKdTreeIndexParams(4); + std::shared_ptr search_params = flannKsTreeSearchParams(32); + + std::shared_ptr index_query = flannKdTreeIndex(points_query, index_params); + std::shared_ptr index_train = flannKdTreeIndex(points_train, index_params); + + // для каждой точки найти total neighbors ближайших соседей + cv::Mat indices_query(n_matches, total_neighbours, CV_32SC1); + cv::Mat distances2_query(n_matches, total_neighbours, CV_32FC1); + cv::Mat indices_train(n_matches, total_neighbours, CV_32SC1); + cv::Mat distances2_train(n_matches, total_neighbours, CV_32FC1); + + index_query->knnSearch(points_query, indices_query, distances2_query, total_neighbours, *search_params); + index_train->knnSearch(points_train, indices_train, distances2_train, total_neighbours, *search_params); + + // оценить радиус поиска для каждой картинки + // NB: radius2_query, radius2_train: квадраты радиуса! + float radius2_query, radius2_train; + { + std::vector max_dists2_query(n_matches); + std::vector max_dists2_train(n_matches); + for (int i = 0; i < n_matches; ++i) { + max_dists2_query[i] = distances2_query.at(i, total_neighbours - 1); + max_dists2_train[i] = distances2_train.at(i, total_neighbours - 1); + } + + int median_pos = n_matches / 2; + std::nth_element(max_dists2_query.begin(), max_dists2_query.begin() + median_pos, max_dists2_query.end()); + std::nth_element(max_dists2_train.begin(), max_dists2_train.begin() + median_pos, max_dists2_train.end()); + + radius2_query = max_dists2_query[median_pos] * radius_limit_scale * radius_limit_scale; + radius2_train = max_dists2_train[median_pos] * radius_limit_scale * radius_limit_scale; + } + + // метч остается, если левое и правое множества первых total_neighbors соседей в радиусах поиска(radius2_query, radius2_train) имеют как минимум consistent_matches общих элементов + // TODO заполнить filtered_matches + filtered_matches.reserve(n_matches); + for(int i=0; i < n_matches; i++ ){ + unordered_set closest_points; + int points_counter =0; + for(int j=0; j(i,j)<= radius2_query){ + closest_points.insert(indices_query.at(i,j)); + } + } + for(int j=0; j(i,j)<= radius2_train) && (closest_points.find(indices_train.at(i,j))) !=closest_points.end()){ + points_counter++; + } + } + if (points_counter>= 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..fe0fb9cf 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,16 @@ 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 dists(query_desc.rows, k, CV_32FC1); + flann_index->knnSearch(query_desc, indices, dists,k, *search_params); + matches.resize(query_desc.rows); + for(int qd_idx = 0; qd_idx< query_desc.rows; qd_idx++){ + std::vector &d_matches = matches[qd_idx]; + d_matches.resize(k); + for (int i =0; i< k; i++){ + cv::DMatch match(qd_idx, indices.at(qd_idx, i),0,std::sqrt(dists.at(qd_idx, i))); + d_matches[i] = match; + } + } } diff --git a/src/phg/sfm/homography.cpp b/src/phg/sfm/homography.cpp index 5cbc780c..22ffbe18 100644 --- a/src/phg/sfm/homography.cpp +++ b/src/phg/sfm/homography.cpp @@ -84,8 +84,17 @@ 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 +177,57 @@ namespace { // * (простое описание для понимания) // * [3] http://ikrisoft.blogspot.com/2015/01/ransac-with-contrario-approach.html -// const int n_matches = points_lhs.size(); -// -// // https://en.wikipedia.org/wiki/Random_sample_consensus#Parameters -// const int n_trials = TODO; -// -// const int n_samples = TODO; -// uint64_t seed = 1; -// const double reprojection_error_threshold_px = 2; -// -// int best_support = 0; -// cv::Mat best_H; -// -// std::vector sample; -// for (int i_trial = 0; i_trial < n_trials; ++i_trial) { -// randomSample(sample, n_matches, n_samples, &seed); -// -// cv::Mat H = estimateHomography4Points(points_lhs[sample[0]], points_lhs[sample[1]], points_lhs[sample[2]], points_lhs[sample[3]], -// points_rhs[sample[0]], points_rhs[sample[1]], points_rhs[sample[2]], points_rhs[sample[3]]); -// -// int support = 0; -// for (int i_point = 0; i_point < n_matches; ++i_point) { -// try { -// cv::Point2d proj = phg::transformPoint(points_lhs[i_point], H); -// if (cv::norm(proj - cv::Point2d(points_rhs[i_point])) < reprojection_error_threshold_px) { -// ++support; -// } -// } catch (const std::exception &e) -// { -// std::cerr << e.what() << std::endl; -// } -// } -// -// if (support > best_support) { -// best_support = support; -// best_H = H; -// -// std::cout << "estimateHomographyRANSAC : support: " << best_support << "/" << n_matches << std::endl; -// -// if (best_support == n_matches) { -// break; -// } -// } -// } -// -// std::cout << "estimateHomographyRANSAC : best support: " << best_support << "/" << n_matches << std::endl; -// -// if (best_support == 0) { -// throw std::runtime_error("estimateHomographyRANSAC : failed to estimate homography"); -// } -// -// return best_H; + const int n_matches = points_lhs.size(); + + // https://en.wikipedia.org/wiki/Random_sample_consensus#Parameters + const int n_trials = 2000; + + 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 +247,15 @@ 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 Td; + T.convertTo(Td, CV_64F); + + const double x = pt.x, y = pt.y; + const double z = Td.at(2, 0) * x + Td.at(2, 1) * y + Td.at(2, 2); + if (std::abs(z) < 1e-8) + throw std::runtime_error("transformPoint: unacceptably small scale"); + + return {(Td.at(0, 0) * x + Td.at(0, 1) * y + Td.at(0, 2)) / z, (Td.at(1, 0) * x + Td.at(1, 1) * y + Td.at(1, 2)) / z}; } 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..2697ef8b 100644 --- a/src/phg/sfm/panorama_stitcher.cpp +++ b/src/phg/sfm/panorama_stitcher.cpp @@ -4,6 +4,8 @@ #include #include +#include + /* * imgs - список картинок * parent - список индексов, каждый индекс указывает, к какой картинке должна быть приклеена текущая картинка @@ -20,11 +22,32 @@ 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"); +{ + int root; + std::vector> children(n_images); + for (int i = 0; i < n_images; ++i) { + if (parent[i] != -1) { + children[parent[i]].push_back(i); + }else{ + root = i; + } } + + Hs[root] = cv::Mat::eye(3, 3, CV_64F); + std::queue q; + q.push(root); + + while (!q.empty()) { + int current = q.front(); + q.pop(); + + for (int child : children[current]) { + cv::Mat H_child_to_current = homography_builder(imgs[child], imgs[current]); + Hs[child] = H_child_to_current * Hs[current]; + q.push(child); + } + } +} bbox2 bbox; for (int i = 0; i < n_images; ++i) {