From 1e26593c4b4151a99adc91496936f44ca023e86c Mon Sep 17 00:00:00 2001 From: Ivan Osokin Date: Tue, 24 Feb 2026 00:01:14 +0300 Subject: [PATCH 1/3] add task01 --- .gitignore | 1 + src/phg/sift/sift.cpp | 189 +++++++++++++++++++++++++----------------- tests/test_sift.cpp | 2 +- 3 files changed, 115 insertions(+), 77 deletions(-) diff --git a/.gitignore b/.gitignore index 54b965f..8e6a2e6 100755 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,4 @@ .idea build cmake-build* +.cache \ No newline at end of file diff --git a/src/phg/sift/sift.cpp b/src/phg/sift/sift.cpp index 7204771..e45a350 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,25 @@ 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::Mat downscaled; + cv::resize(oct.layers[s], downscaled, cv::Size(), 0.5, 0.5, cv::INTER_NEAREST); + base = std::move(downscaled); // можно использовать и downsample2x_avg(oct.layers[s]), это позволяет потом заапскейлить слои обратно до оригинального разрешения без сдвига // но потребуется везде изменить формулу для пересчета ключевых точек: pt_upscaled = (pt_downscaled + 0.5) * 2^o - 0.5 @@ -139,6 +151,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 +265,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 +314,17 @@ std::vector phg::findScaleSpaceExtrema(const std::vector (r + 1) * (r + 1) * det) { + break; + } } // скейлим координаты точек обратно до родных размеров картинки @@ -379,39 +416,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 +478,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 +611,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 +646,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 +658,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_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 From c04abb8825c3eb55b4ba264abc3689067c4edbaf Mon Sep 17 00:00:00 2001 From: Ivan Osokin Date: Thu, 26 Feb 2026 22:06:45 +0300 Subject: [PATCH 2/3] small fix --- src/phg/sift/sift.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/phg/sift/sift.cpp b/src/phg/sift/sift.cpp index e45a350..9f5fcb9 100755 --- a/src/phg/sift/sift.cpp +++ b/src/phg/sift/sift.cpp @@ -128,9 +128,7 @@ std::vector phg::buildOctaves(const cv::Mat& img, const phg:: if (o + 1 < n_octaves) { // используется в opencv, формула для пересчета ключевых точек: pt_upscaled = 2^o * pt_downscaled // TODO cv::resize(даунскейлим текущий слой в два раза, без интерполяции, просто сабсепмлинг); - cv::Mat downscaled; - cv::resize(oct.layers[s], downscaled, cv::Size(), 0.5, 0.5, cv::INTER_NEAREST); - base = std::move(downscaled); + 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 From 5974f21aeba1bba808bd94afca5edf13f87027a5 Mon Sep 17 00:00:00 2001 From: Ivan Osokin Date: Thu, 26 Feb 2026 23:05:03 +0300 Subject: [PATCH 3/3] bruteforec --- src/phg/sift/sift.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/phg/sift/sift.h b/src/phg/sift/sift.h index 6e65f6d..e19f126 100755 --- a/src/phg/sift/sift.h +++ b/src/phg/sift/sift.h @@ -11,7 +11,7 @@ struct SIFTParams { double edge_threshold = 10; double sigma = 1.6; - double orient_peak_ratio = 0.8; + double orient_peak_ratio = 0.8482222; int orient_nbins = 36; bool upscale_first = true;