From 6f3317ea2b76d9ba6af1e29afc9946f85670517c Mon Sep 17 00:00:00 2001 From: IrDIE <110756720+IrDIE@users.noreply.github.com> Date: Wed, 25 Feb 2026 09:59:20 +0300 Subject: [PATCH] add first solution --- src/phg/sift/sift.cpp | 324 +++++++++++++++++++++++++++++------------- src/phg/sift/sift.h | 6 +- tests/test_sift.cpp | 19 ++- 3 files changed, 240 insertions(+), 109 deletions(-) diff --git a/src/phg/sift/sift.cpp b/src/phg/sift/sift.cpp index 7204771..bc1169a 100755 --- a/src/phg/sift/sift.cpp +++ b/src/phg/sift/sift.cpp @@ -16,6 +16,7 @@ // 3) https://github.com/opencv/opencv/blob/1834eed8098aa2c595f4d1099eeaa0992ce8b321/modules/features2d/src/sift.dispatch.cpp // 4) https://github.com/opencv/opencv/blob/1834eed8098aa2c595f4d1099eeaa0992ce8b321/modules/features2d/src/sift.simd.hpp +// My SIFT namespace { cv::Mat upsample2x(const cv::Mat& src) @@ -88,7 +89,9 @@ 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 пикселей, там уже не имеет смысла что-то детектировать @@ -108,16 +111,23 @@ 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 k = std::pow(2.0, 1.0 / s); + double sigma_layer = sigma0 * std::pow(k, i); + double sigma_prev = sigma0 * std::pow(k, i - 1); + double sigma_adjusted = std::sqrt(sigma_layer * sigma_layer - sigma_prev * sigma_prev); // Incremental sigma + + oct.layers[i] = base.clone(); + cv::GaussianBlur(oct.layers[i - 1], oct.layers[i], cv::Size(), sigma_adjusted, sigma_adjusted); } // подготавливаем базовый слой для следующей октавы if (o + 1 < n_octaves) { + + cv::Mat base_; + cv::resize(oct.layers[n_layers-1], base_, cv::Size(), 0.5, 0.5, cv::INTER_NEAREST); + base = base_; // используется в opencv, формула для пересчета ключевых точек: pt_upscaled = 2^o * pt_downscaled - // TODO cv::resize(даунскейлим текущий слой в два раза, без интерполяции, просто сабсепмлинг); + // T1ODO cv::resize(даунскейлим текущий слой в два раза, без интерполяции, просто сабсепмлинг); // можно использовать и downsample2x_avg(oct.layers[s]), это позволяет потом заапскейлить слои обратно до оригинального разрешения без сдвига // но потребуется везде изменить формулу для пересчета ключевых точек: pt_upscaled = (pt_downscaled + 0.5) * 2^o - 0.5 @@ -131,14 +141,33 @@ std::vector phg::buildOctaves(const cv::Mat& img, const phg:: } std::vector phg::buildDoG(const std::vector& octaves, const phg::SIFTParams& p, int verbose_level) -{ +{ // difference of Gaissians 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); + for (size_t i = 0; i+1 < octave.layers.size(); ++i) { + cv::Mat& dog_layer = dog[o].layers[i]; + const cv::Mat& layer1 = octave.layers[i]; + const cv::Mat& layer2 = octave.layers[i + 1]; + + dog_layer.create(layer1.size(), CV_32F); + + for (int y = 0; y < layer1.rows; y++) { + const float* r1 = layer1.ptr(y); + const float* r2 = layer2.ptr(y); + float* rd = dog_layer.ptr(y); + for (int x = 0; x < layer1.cols; x++) { + rd[x] = r2[x] - r1[x]; + } + } - // TODO каждый слой дога это разница n+1 и n-й гауссианы + } + + + // T1ODO каждый слой дога это разница n+1 и n-й гауссианы } return dog; @@ -193,7 +222,8 @@ std::vector phg::findScaleSpaceExtrema(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, dxy, dyy, dxs, dys, dss; + 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; - cv::Matx33f H(dxx, dxy, dxs, dxy, dyy, dys, dxs, dys, dss); + 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); cv::Vec3f g(dx, dy, ds); // в нашей точке производная (градиент) еще не равна нулю (т.к. еще мы скорее всего не точно в оптимуме) - // хотим найти такой offset, где ноль производной. в предположении что оптимизируемая функция это парабола, ищем корни ее производной, линейной функции + // хотим найти такой offset, где ноль производной. в предположении что оптимизируемая функция это парабола, + // ищем корни ее производной, линейной функции: // grad(x + offset) = grad(x) + grad'(x) * offset = grad(x) + hessian(x) * offset = 0 // hessian(x) * offset = -grad(x) // линейная система. можно решить специализированным решателем либо просто найти обратную матрицу гессиана и домножить на минус градиент // offset = -hessian(x)^-1 * grad(x) - cv::Vec3f offset; if (!cv::solve(H, -g, offset, cv::DECOMP_LU)) break; @@ -272,22 +326,31 @@ std::vector phg::findScaleSpaceExtrema(const std::vector edge_threshold) + break; } // скейлим координаты точек обратно до родных размеров картинки @@ -342,6 +405,7 @@ std::vector phg::findScaleSpaceExtrema(const std::vector phg::computeOrientations(const std::vector& kpts, const std::vector& octaves, const phg::SIFTParams& params, int verbose_level) { + // we wanna assign one or more dominant orientations to each keypoint detected in the scale-space const int s = params.n_octave_layers; const double sigma0 = params.sigma; const int n_bins = params.orient_nbins; @@ -366,8 +430,13 @@ std::vector phg::computeOrientations(const 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; // we take a point within +-rad + // Illustration of gradient computation: + // The pixel (px, py) is at the center, and we compute the gradient + // in the x-direction (gx) and y-direction (gy) using neighboring pixels. + // + // gy (vertical gradient) + // ↑ + // | (py-1, px-1) (py-1, px) (py-1, px+1) + // | +------------+------------+ + // | | | | + // | | | | + // | +------------+------------+ + // | (py, px-1) (py, px) (py, px+1) + // | +------------+------------+ + // | | | | + // | | | | + // | +------------+------------+ + // | (py+1, px-1) (py+1, px) (py+1, px+1) + // | | | | + // | | | | + // | +------------+------------+ + // | + // +--------------------> gx (horizontal gradient) + // + // gx = img.at(py, px+1) - img.at(py, px-1) + // gy = img.at(py+1, px) - img.at(py-1, px) + // gx measures the horizontal change in intensity (difference between right and left neighbors). + // gy measures the vertical change in intensity (difference between top and bottom neighbors). + // Together, gx and gy form the gradient vector at (px, py), which points in the direction of the steepest intensity change. + + // градиент + 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; + if (!params.enable_orientation_gaussian_weighting) { + weight = 1.f; + } else { + weight = std::exp(-(dx*dx + dy*dy) / (2.f * sigma_win * sigma_win)); + } + + // голосуем в гистограмме направлений. находим два ближайших бина и гладко распределяем голос между ними + // в таком случае, голос попавший близко к границе между бинами, проголосует поровну за оба бина + // float bins_closest = angle_deg / (360./n_bins); + float bin = angle_deg / (360./n_bins); + if (bin >= n_bins) bin -= n_bins; + int bin0 = (int) bin; + int bin1 = (bin0 + 1) % n_bins; + + float frac = bin - bin0; + if (!params.enable_orientation_bin_interpolation) { + frac = 0.f; + } + + histogram[bin0] += weight * mag * (1. - frac); + histogram[bin1] += weight * mag * (frac); } } + + // немного сгладим гистограмму: сделаем несколько проходов box-blur (повторный box blur приближает гауссово размытие) for (int iter = 0; iter < 6; iter++) { float first = histogram[0]; @@ -428,7 +530,11 @@ std::vector phg::computeOrientations(const std::vector phg::computeOrientations(const std::vector a = (left + right - 2 * center) / 2 // f(1) - f(-1) = 2b -> b = (right - left) / 2 + 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 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 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); } } } @@ -560,6 +669,22 @@ std::pair> phg::computeDescriptors(const std: float mag = std::sqrt(gx * gx + gy * gy); float angle = std::atan2(gy, gx); + /* I had a question about why we need gradient both for computeOrientation and computeDescriptors: + + Why Can't We Reuse the Gradients and Angles? + Different Neighborhoods: + In computeOrientations, the gradients are computed over a circular region around the keypoint to determine the dominant orientation. + In computeDescriptors, the gradients are computed over a larger rectangular region divided into spatial bins to build the descriptor. + + Different Uses: + In computeOrientations, the gradients are used to build an orientation histogram. + In computeDescriptors, the gradients are used to populate the spatial and orientation bins of the descriptor. + + Rotation Adjustment: + In computeDescriptors, the gradient angles are adjusted relative to the keypoint's dominant orientation. This adjustment is not needed in computeOrientations. + Performance Considerations: + + */ // инвариантность к повороту: повернем направление градиента на угол ключевой точки float angle_invariant = angle - kp_angle_rad; @@ -574,11 +699,12 @@ 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(-(dx*dx + dy*dy) / (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 +735,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] += wo * wx * wy * weighted_mag; } } } @@ -620,9 +746,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[idx] += weighted_mag; } } } @@ -730,7 +855,8 @@ void phg::SIFT::detectAndCompute(const cv::Mat& img, const cv::Mat& mask, std::v std::tie(desc, kpts) = computeDescriptors(kpts, octaves, p, verbose_level); - // TODO всегда ли мы получаем ровно столько точек сколько запросили в параметре nfeatures? в каких случаях это не так и в какую сторону? + // TODO всегда ли мы получаем ровно столько точек сколько запросили в параметре nfeatures? + // в каких случаях это не так и в какую сторону? // как подкрутить алгоритм, чтобы всегда выдавать ровно запрошенное количество точек (когда это в принципе возможно) но не сильно просесть в производительности? } diff --git a/src/phg/sift/sift.h b/src/phg/sift/sift.h index 6e65f6d..0dc3647 100755 --- a/src/phg/sift/sift.h +++ b/src/phg/sift/sift.h @@ -3,7 +3,7 @@ #include namespace phg { - +// rotation 30 :: [ORB_OCV] average angle difference between matched points: 25.313 degrees struct SIFTParams { int nfeatures = 0; int n_octave_layers = 3; @@ -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.5; int orient_nbins = 36; bool upscale_first = true; @@ -30,7 +30,7 @@ class SIFT { std::vector layers; }; - explicit SIFT(const SIFTParams& p = SIFTParams(), int verbose_level = 0, const std::string& debug_folder = "") + explicit SIFT(const SIFTParams& p = SIFTParams(), int verbose_level = 0, const std::string& debug_folder = "/mnt/c/Users/irady/GitHub/ITMO/2026-фотграмметрия/homeworks_wsl/task01/PhotogrammetryTasks2026/data/debug/test_sift/debug/") : p(p) , verbose_level(verbose_level) , debug_folder(debug_folder) diff --git a/tests/test_sift.cpp b/tests/test_sift.cpp index cf3bd7d..8d8b1b5 100755 --- a/tests/test_sift.cpp +++ b/tests/test_sift.cpp @@ -25,13 +25,15 @@ #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 + +// [ORB_OCV] average angle difference between matched points: 25.313 degrees +// [SIFTOCV] average angle difference between matched points: 9.58379 degrees with SIFT_ORI_PEAK_RATIO 0.8 +// [SIFTOCV] average angle difference between matched points: 24.1679 degrees with SIFT_ORI_PEAK_RATIO 0.999 struct MatchingPairData { size_t npoints1, npoints2, nmatches; }; @@ -149,7 +151,8 @@ void evaluateDetection(const cv::Mat& M, double minRecall, cv::Mat img0 = cv::Ma size_t width = img0.cols; size_t height = img0.rows; cv::Mat transformedImage; - cv::warpAffine(img0, transformedImage, M, cv::Size(width, height)); // строим img1 - преобразованная исходная картинка в соответствии с закодированным в матрицу M искажением пространства + cv::warpAffine(img0, transformedImage, M, cv::Size(width, height)); // строим img1 - преобразованная исходная картинка в соответствии + // с закодированным в матрицу M искажением пространства cv::Mat noise(cv::Size(width, height), CV_8UC3); cv::setRNGSeed(125125); // фиксируем рандом для детерминизма (чтобы результат воспроизводился из раза в раз) cv::randn(noise, cv::Scalar::all(0), cv::Scalar::all(GAUSSIAN_NOISE_STDDEV)); @@ -160,6 +163,7 @@ void evaluateDetection(const cv::Mat& M, double minRecall, cv::Mat img0 = cv::Ma for (int method = 0; method < 3; ++method) { // тестируем три метода: OpenCV ORB, OpenCV SIFT, ваш SIFT std::vector kps0; std::vector kps1; + // size ~ diameter of the meaningful neighborhood around that point cv::Mat desc0; cv::Mat desc1; @@ -178,7 +182,7 @@ void evaluateDetection(const cv::Mat& M, double minRecall, cv::Mat img0 = cv::Ma detector->detect(img0, kps0); // детектируем ключевые точки на исходной картинке detector->detect(img1, kps1); // детектируем ключевые точки на преобразованной картинке - detector->compute(img0, kps0, desc0); + detector->compute(img0, kps0, desc0); // compute descriptors for kpts detector->compute(img1, kps1, desc1); } else if (method == 1) { method_name = "SIFTOCV"; @@ -304,7 +308,8 @@ void evaluateDetection(const cv::Mat& M, double minRecall, cv::Mat img0 = cv::Ma std::cout << log_prefix << "average size ratio between matched points: " << (size_ratio_sum / n_matched) << std::endl; if (angle_diff_sum != 0.0) { std::cout << log_prefix << "average angle difference between matched points: " << (angle_diff_sum / n_matched) << " degrees" << std::endl; - // TODO почему SIFT менее точно угадывает средний угол отклонения? изменяется ли ситуация если выкрутить параметр ORIENTATION_VOTES_PEAK_RATIO=0.999? почему? + // TODO почему SIFT менее точно угадывает средний угол отклонения? + // изменяется ли ситуация если выкрутить параметр ORIENTATION_VOTES_PEAK_RATIO=0.999? почему? } if (desc_dist_sum != 0.0 && desc_rand_dist_sum != 0.0) { std::cout << log_prefix << "average descriptor distance between matched points: " << (desc_dist_sum / n_matched) << " (random distance: " << (desc_rand_dist_sum / n_matched)