Skip to content

Task01 Иван Осокин ИТМО #7

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
ivan0sokin wants to merge 3 commits into from
+114 −78
Closed

Task01 Иван Осокин ИТМО #7

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
1e26593
add task01
Feb 23, 2026
c04abb8
small fix
Feb 26, 2026
5974f21
bruteforec
Feb 26, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,3 +2,4 @@
.idea
build
cmake-build*
.cache
187 changes: 111 additions & 76 deletions src/phg/sift/sift.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,8 +1,12 @@
#include "sift.h"
#include "libutils/rasserts.h"

#include <cmath>
#include <exception>
#include <iostream>
#include <numeric>
#include <opencv2/core/mat.hpp>
#include <opencv2/core/types.hpp>
#include <opencv2/features2d.hpp>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/imgproc.hpp>
Expand Down Expand Up @@ -79,6 +83,7 @@ cv::Mat phg::toGray32F(const cv::Mat& img)
std::vector<phg::SIFT::Octave> 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 не соглашается со "взятое с потолка":
Expand Down Expand Up @@ -107,18 +112,23 @@ std::vector<phg::SIFT::Octave> 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::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

Expand All @@ -139,6 +149,9 @@ std::vector<phg::SIFT::Octave> phg::buildDoG(const std::vector<phg::SIFT::Octave
dog[o].layers.resize(octave.layers.size() - 1);

// TODO каждый слой дога это разница n+1 и n-й гауссианы
for (int i = 0; i < dog[o].layers.size(); i++) {
dog[o].layers[i] = octave.layers[i + 1] - octave.layers[i];
}
}

return dog;
Expand Down Expand Up @@ -206,19 +219,31 @@ std::vector<cv::KeyPoint> phg::findScaleSpaceExtrema(const std::vector<phg::SIFT
if (v <= val)
is_min = false;
};

auto checkNeighbors = [&](const float* prev, const float* curr, const float* next, const std::vector<int> 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;

Expand All @@ -238,13 +263,12 @@ std::vector<cv::KeyPoint> phg::findScaleSpaceExtrema(const std::vector<phg::SIFT

// гессиан
float dxx, dxy, dyy, dxs, dys, dss;
// float dxx = cL.at<float>(yi, xi + 1) + cL.at<float>(yi, xi - 1) - 2.f * resp_center;
// float dyy = TODO;
// float dss = TODO;
//
// float dxy = (cL.at<float>(yi + 1, xi + 1) - cL.at<float>(yi + 1, xi - 1) - cL.at<float>(yi - 1, xi + 1) + cL.at<float>(yi - 1, xi - 1)) * 0.25f;
// float dxs = TODO;
// float dys = TODO;
dxx = cL.at<float>(yi, xi + 1) + cL.at<float>(yi, xi - 1) - 2.f * resp_center;
dyy = cL.at<float>(yi + 1, xi) + cL.at<float>(yi - 1, xi) - 2.f * resp_center;
dss = nL.at<float>(yi, xi) + pL.at<float>(yi, xi) - 2.f * resp_center;
dxy = (cL.at<float>(yi + 1, xi + 1) - cL.at<float>(yi + 1, xi - 1) - cL.at<float>(yi - 1, xi + 1) + cL.at<float>(yi - 1, xi - 1)) * 0.25f;
dxs = (nL.at<float>(yi, xi + 1) - nL.at<float>(yi, xi - 1) - pL.at<float>(yi, xi + 1) + pL.at<float>(yi, xi - 1)) * 0.25f;
dys = (nL.at<float>(yi + 1, xi) - nL.at<float>(yi - 1, xi) - pL.at<float>(yi + 1, xi) + pL.at<float>(yi - 1, xi)) * 0.25f;

cv::Matx33f H(dxx, dxy, dxs, dxy, dyy, dys, dxs, dys, dss);

Expand Down Expand Up @@ -288,6 +312,17 @@ std::vector<cv::KeyPoint> phg::findScaleSpaceExtrema(const std::vector<phg::SIFT
// float r = edge_threshold;
// if (TODO)
// break;

float trace = dxx + dyy;
float det = dxx * dyy - dxy * dxy;
if (det <= 0) {
break;
}

float r = edge_threshold;
if (trace * trace * r > (r + 1) * (r + 1) * det) {
break;
}
}

// скейлим координаты точек обратно до родных размеров картинки
Expand Down Expand Up @@ -379,39 +414,39 @@ std::vector<cv::KeyPoint> phg::computeOrientations(const std::vector<cv::KeyPoin

for (int dy = -radius; dy <= radius; dy++) {
for (int dx = -radius; dx <= radius; dx++) {
// int px = xi + dx;
// int py = yi + dy;
//
// // градиент
// float gx = img.at<float>(py, px + 1) - img.at<float>(py, px - 1);
// float gy = img.at<float>(py + 1, px) - img.at<float>(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<float>(py, px + 1) - img.at<float>(py, px - 1);
float gy = img.at<float>(py + 1, px) - img.at<float>(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;
}
}

Expand Down Expand Up @@ -441,29 +476,29 @@ std::vector<cv::KeyPoint> phg::computeOrientations(const std::vector<cv::KeyPoin
float center = histogram[i];
float right = histogram[next];

// хотим найти угол дескриптора точнее = зафитить параболу по трем точкам (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 = 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);
}
}
}
Expand Down Expand Up @@ -574,11 +609,11 @@ std::pair<cv::Mat, std::vector<cv::KeyPoint>> 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 по пространственным бинам и по бинам гистограммок трилинейной интерполяцией
Expand Down Expand Up @@ -609,8 +644,8 @@ std::pair<cv::Mat, std::vector<cv::KeyPoint>> 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;
}
}
}
Expand All @@ -621,8 +656,8 @@ std::pair<cv::Mat, std::vector<cv::KeyPoint>> 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;
}
}
}
Expand Down
2 changes: 1 addition & 1 deletion src/phg/sift/sift.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down
2 changes: 1 addition & 1 deletion tests/test_sift.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
Loading