From a1ce9c07f0310fee0b1cfc7a94417b79cb110dc0 Mon Sep 17 00:00:00 2001 From: NovokshanovE Date: Sat, 28 Feb 2026 22:05:44 +0300 Subject: [PATCH 1/3] Try #1 --- src/phg/sift/sift.cpp | 175 ++++++++++++++++++++++++++++-------------- tests/test_sift.cpp | 2 +- 2 files changed, 118 insertions(+), 59 deletions(-) diff --git a/src/phg/sift/sift.cpp b/src/phg/sift/sift.cpp index 72047711..b05582f4 100755 --- a/src/phg/sift/sift.cpp +++ b/src/phg/sift/sift.cpp @@ -111,13 +111,16 @@ std::vector phg::buildOctaves(const cv::Mat& img, const phg:: // 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); + 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(даунскейлим текущий слой в два раза, без интерполяции, просто сабсепмлинг); + cv::resize(oct.layers[s], base, cv::Size(oct.layers[s].cols / 2, oct.layers[s].rows / 2), 0, 0, cv::INTER_NEAREST); // можно использовать и downsample2x_avg(oct.layers[s]), это позволяет потом заапскейлить слои обратно до оригинального разрешения без сдвига // но потребуется везде изменить формулу для пересчета ключевых точек: pt_upscaled = (pt_downscaled + 0.5) * 2^o - 0.5 @@ -139,6 +142,9 @@ std::vector phg::buildDoG(const std::vector phg::findScaleSpaceExtrema(const std::vector phg::findScaleSpaceExtrema(const std::vector(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; // float dxx = cL.at(yi, xi + 1) + cL.at(yi, xi - 1) - 2.f * resp_center; // float dyy = TODO; // float dss = TODO; @@ -273,6 +311,21 @@ std::vector phg::findScaleSpaceExtrema(const std::vector= (r + 1.f) * (r + 1.f) * det) + break; // float trace = //TODO ; // = lambda1 + lambda2 // float det = // TODO ; // = lambda1 * lambda2 // if (det <= 0) @@ -379,39 +432,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::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] += mag * weight * (1.f - frac); + histogram[bin1] += mag * weight * frac; } } @@ -448,22 +503,26 @@ std::vector phg::computeOrientations(const std::vector a = (left + right - 2 * center) / 2 - // f(1) - f(-1) = 2b -> b = (right - left) / 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 denom = left + right - 2.f * center; + float offset = (std::abs(denom) > 1e-7f) ? (0.5f * (left - right) / denom) : 0.f; + if (!params.enable_orientation_subpixel_localization) { + offset = 0.f; + } + offset = std::max(-1.f, std::min(1.f, offset)); + + 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 +633,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 +668,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 +680,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 cf3bd7df..7433f9bb 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 2a63617f930e6e949752e2123f3ee45aba39f832 Mon Sep 17 00:00:00 2001 From: NovokshanovE Date: Tue, 3 Mar 2026 11:03:41 +0300 Subject: [PATCH 2/3] =?UTF-8?q?=D0=94=D0=BE=D0=B1=D0=B0=D0=B2=D0=B8=D0=BB?= =?UTF-8?q?=20OpenMP?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../googletest/src/gtest-death-test.cc | 1 + libs/utils/libutils/thread_mutex.h | 1 + src/phg/sift/sift.cpp | 67 +++++++++++++++---- 3 files changed, 56 insertions(+), 13 deletions(-) diff --git a/libs/3rdparty/libgtest/googletest/src/gtest-death-test.cc b/libs/3rdparty/libgtest/googletest/src/gtest-death-test.cc index 43cbd78a..36626fcd 100755 --- a/libs/3rdparty/libgtest/googletest/src/gtest-death-test.cc +++ b/libs/3rdparty/libgtest/googletest/src/gtest-death-test.cc @@ -45,6 +45,7 @@ # include # include +# include # include # if GTEST_OS_LINUX diff --git a/libs/utils/libutils/thread_mutex.h b/libs/utils/libutils/thread_mutex.h index 25f4e841..b3ebf74c 100755 --- a/libs/utils/libutils/thread_mutex.h +++ b/libs/utils/libutils/thread_mutex.h @@ -91,6 +91,7 @@ class TryLock { bool acquire () { _locked = _mutex.tryLock(); + return _locked; } void release () diff --git a/src/phg/sift/sift.cpp b/src/phg/sift/sift.cpp index b05582f4..91cc960a 100755 --- a/src/phg/sift/sift.cpp +++ b/src/phg/sift/sift.cpp @@ -107,6 +107,7 @@ std::vector phg::buildOctaves(const cv::Mat& img, const phg:: // для простоты в каждой октаве будем каждый раз блюрить базовую картинку с полной сигмой // можно подумать, как сделать эффективнее - для построения n+1 слоя доблюревать уже поблюренный n-ый слой, так чтобы в итоге получилась такая же сигма // это будет немного быстрее, тк нужно более маленькое ядро свертки на каждый шаг + #pragma omp parallel for if (n_layers > 3) for (int i = 1; i < n_layers; i++) { // TODO double sigma_layer = sigma0 * корень из двух нужной степени, чтобы при i==s получали удвоение базового блюра; // // вычтем sigma0 чтобы размыть ровно до нужной суммарной сигмы @@ -137,6 +138,7 @@ std::vector phg::buildDoG(const std::vector dog(octaves.size()); + #pragma omp parallel for if (octaves.size() > 1) for (size_t o = 0; o < octaves.size(); o++) { const phg::SIFT::Octave& octave = octaves[o]; dog[o].layers.resize(octave.layers.size() - 1); @@ -166,14 +168,16 @@ std::vector phg::findScaleSpaceExtrema(const std::vector keypoints; + std::vector> keypoints_by_octave(dog.size()); + #pragma omp parallel for if (dog.size() > 1) 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(); rassert(n_dog_layers == s + 2, 2138971238612312); + std::vector& keypoints = keypoints_by_octave[o]; // итерируемся по внутренним слоям пирамиды, у нас всегда есть предыдущий и следующий сосед for (int layer = 1; layer <= s; layer++) { @@ -384,7 +388,16 @@ std::vector phg::findScaleSpaceExtrema(const std::vector keypoints; + size_t total_keypoints = 0; + for (const auto& octave_keypoints : keypoints_by_octave) + total_keypoints += octave_keypoints.size(); + keypoints.reserve(total_keypoints); + for (const auto& octave_keypoints : keypoints_by_octave) { + keypoints.insert(keypoints.end(), octave_keypoints.begin(), octave_keypoints.end()); } if (verbose_level) @@ -400,13 +413,13 @@ std::vector phg::computeOrientations(const std::vector histogram(n_bins); - - std::vector oriented_kpts; + std::vector> oriented_by_kpt(kpts.size()); const int first_octave = params.upscale_first ? -1 : 0; - for (const cv::KeyPoint& kp : kpts) { + #pragma omp parallel for if (kpts.size() > 32) + for (int kp_idx = 0; kp_idx < (int)kpts.size(); kp_idx++) { + const cv::KeyPoint& kp = kpts[kp_idx]; int layer = kp.class_id; int real_octave = kp.octave; int o = real_octave - first_octave; // индекс в массиве octaves @@ -428,6 +441,7 @@ std::vector phg::computeOrientations(const std::vector= img.cols - 1 || yi - radius <= 0 || yi + radius >= img.rows - 1) continue; + std::vector histogram(n_bins); histogram.assign(n_bins, 0.0); for (int dy = -radius; dy <= radius; dy++) { @@ -522,11 +536,20 @@ std::vector phg::computeOrientations(const std::vector oriented_kpts; + size_t total_oriented = 0; + for (const auto& kp_orientations : oriented_by_kpt) + total_oriented += kp_orientations.size(); + oriented_kpts.reserve(total_oriented); + for (const auto& kp_orientations : oriented_by_kpt) { + oriented_kpts.insert(oriented_kpts.end(), kp_orientations.begin(), kp_orientations.end()); + } + if (verbose_level) std::cout << "orientations: " << kpts.size() << " -> " << oriented_kpts.size() << " keypoints" << std::endl; @@ -555,7 +578,12 @@ std::pair> phg::computeDescriptors(const std: const int first_octave = params.upscale_first ? -1 : 0; - for (const cv::KeyPoint& kp : kpts) { + std::vector> descriptor_rows(kpts.size()); + std::vector descriptor_valid(kpts.size(), 0); + + #pragma omp parallel for if (kpts.size() > 32) + for (int kp_idx = 0; kp_idx < (int)kpts.size(); kp_idx++) { + const cv::KeyPoint& kp = kpts[kp_idx]; int layer = kp.class_id; int real_octave = kp.octave; int o = real_octave - first_octave; // индекс в массиве octaves @@ -707,13 +735,26 @@ std::pair> phg::computeDescriptors(const std: for (float& v : desc) v /= norm; - if (descriptors.empty()) { - descriptors.create(0, n_dims, CV_32F); - } + descriptor_rows[kp_idx] = std::move(desc); + descriptor_valid[kp_idx] = 1; + } + + size_t valid_count = 0; + for (unsigned char is_valid : descriptor_valid) + valid_count += is_valid; + + if (valid_count > 0) { + descriptors.create(0, n_dims, CV_32F); + valid_kpts.reserve(valid_count); + } + + for (int kp_idx = 0; kp_idx < (int)kpts.size(); kp_idx++) { + if (!descriptor_valid[kp_idx]) + continue; - cv::Mat row(1, n_dims, CV_32F, desc.data()); + cv::Mat row(1, n_dims, CV_32F, descriptor_rows[kp_idx].data()); descriptors.push_back(row.clone()); - valid_kpts.push_back(kp); + valid_kpts.push_back(kpts[kp_idx]); } if (verbose_level) From b9b506ce3a1784047173aedc666e2f221df15fb3 Mon Sep 17 00:00:00 2001 From: NovokshanovE Date: Tue, 3 Mar 2026 11:58:29 +0300 Subject: [PATCH 3/3] =?UTF-8?q?=D0=A3=D0=B1=D1=80=D0=B0=D0=BB=20openMP?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/phg/sift/sift.cpp | 68 +++++++++---------------------------------- 1 file changed, 13 insertions(+), 55 deletions(-) diff --git a/src/phg/sift/sift.cpp b/src/phg/sift/sift.cpp index 91cc960a..1eeda8ae 100755 --- a/src/phg/sift/sift.cpp +++ b/src/phg/sift/sift.cpp @@ -107,7 +107,6 @@ std::vector phg::buildOctaves(const cv::Mat& img, const phg:: // для простоты в каждой октаве будем каждый раз блюрить базовую картинку с полной сигмой // можно подумать, как сделать эффективнее - для построения n+1 слоя доблюревать уже поблюренный n-ый слой, так чтобы в итоге получилась такая же сигма // это будет немного быстрее, тк нужно более маленькое ядро свертки на каждый шаг - #pragma omp parallel for if (n_layers > 3) for (int i = 1; i < n_layers; i++) { // TODO double sigma_layer = sigma0 * корень из двух нужной степени, чтобы при i==s получали удвоение базового блюра; // // вычтем sigma0 чтобы размыть ровно до нужной суммарной сигмы @@ -138,7 +137,6 @@ std::vector phg::buildDoG(const std::vector dog(octaves.size()); - #pragma omp parallel for if (octaves.size() > 1) for (size_t o = 0; o < octaves.size(); o++) { const phg::SIFT::Octave& octave = octaves[o]; dog[o].layers.resize(octave.layers.size() - 1); @@ -168,16 +166,13 @@ std::vector phg::findScaleSpaceExtrema(const std::vector> keypoints_by_octave(dog.size()); - - #pragma omp parallel for if (dog.size() > 1) + std::vector keypoints; 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(); rassert(n_dog_layers == s + 2, 2138971238612312); - std::vector& keypoints = keypoints_by_octave[o]; // итерируемся по внутренним слоям пирамиды, у нас всегда есть предыдущий и следующий сосед for (int layer = 1; layer <= s; layer++) { @@ -388,16 +383,7 @@ std::vector phg::findScaleSpaceExtrema(const std::vector keypoints; - size_t total_keypoints = 0; - for (const auto& octave_keypoints : keypoints_by_octave) - total_keypoints += octave_keypoints.size(); - keypoints.reserve(total_keypoints); - for (const auto& octave_keypoints : keypoints_by_octave) { - keypoints.insert(keypoints.end(), octave_keypoints.begin(), octave_keypoints.end()); + std::cout << "octave " << o << ": " << keypoints.size() << " keypoints so far" << std::endl; } if (verbose_level) @@ -413,13 +399,13 @@ std::vector phg::computeOrientations(const std::vector> oriented_by_kpt(kpts.size()); + std::vector histogram(n_bins); + + std::vector oriented_kpts; const int first_octave = params.upscale_first ? -1 : 0; - #pragma omp parallel for if (kpts.size() > 32) - for (int kp_idx = 0; kp_idx < (int)kpts.size(); kp_idx++) { - const cv::KeyPoint& kp = kpts[kp_idx]; + for (const cv::KeyPoint& kp : kpts) { int layer = kp.class_id; int real_octave = kp.octave; int o = real_octave - first_octave; // индекс в массиве octaves @@ -441,7 +427,6 @@ std::vector phg::computeOrientations(const std::vector= img.cols - 1 || yi - radius <= 0 || yi + radius >= img.rows - 1) continue; - std::vector histogram(n_bins); histogram.assign(n_bins, 0.0); for (int dy = -radius; dy <= radius; dy++) { @@ -536,20 +521,11 @@ std::vector phg::computeOrientations(const std::vector oriented_kpts; - size_t total_oriented = 0; - for (const auto& kp_orientations : oriented_by_kpt) - total_oriented += kp_orientations.size(); - oriented_kpts.reserve(total_oriented); - for (const auto& kp_orientations : oriented_by_kpt) { - oriented_kpts.insert(oriented_kpts.end(), kp_orientations.begin(), kp_orientations.end()); - } - if (verbose_level) std::cout << "orientations: " << kpts.size() << " -> " << oriented_kpts.size() << " keypoints" << std::endl; @@ -578,12 +554,7 @@ std::pair> phg::computeDescriptors(const std: const int first_octave = params.upscale_first ? -1 : 0; - std::vector> descriptor_rows(kpts.size()); - std::vector descriptor_valid(kpts.size(), 0); - - #pragma omp parallel for if (kpts.size() > 32) - for (int kp_idx = 0; kp_idx < (int)kpts.size(); kp_idx++) { - const cv::KeyPoint& kp = kpts[kp_idx]; + for (const cv::KeyPoint& kp : kpts) { int layer = kp.class_id; int real_octave = kp.octave; int o = real_octave - first_octave; // индекс в массиве octaves @@ -735,26 +706,13 @@ std::pair> phg::computeDescriptors(const std: for (float& v : desc) v /= norm; - descriptor_rows[kp_idx] = std::move(desc); - descriptor_valid[kp_idx] = 1; - } - - size_t valid_count = 0; - for (unsigned char is_valid : descriptor_valid) - valid_count += is_valid; - - if (valid_count > 0) { - descriptors.create(0, n_dims, CV_32F); - valid_kpts.reserve(valid_count); - } - - for (int kp_idx = 0; kp_idx < (int)kpts.size(); kp_idx++) { - if (!descriptor_valid[kp_idx]) - continue; + if (descriptors.empty()) { + descriptors.create(0, n_dims, CV_32F); + } - cv::Mat row(1, n_dims, CV_32F, descriptor_rows[kp_idx].data()); + cv::Mat row(1, n_dims, CV_32F, desc.data()); descriptors.push_back(row.clone()); - valid_kpts.push_back(kpts[kp_idx]); + valid_kpts.push_back(kp); } if (verbose_level)