Skip to content

Task01 Вадим Бенкевич НИИЧАВО #10

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

Open
vb002 wants to merge 2 commits into
base: task01
Choose a base branch
from
Open

Task01 Вадим Бенкевич НИИЧАВО #10

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
33e7abf
Drop code
Feb 25, 2026
16779da
Fill the gaps
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
159 changes: 103 additions & 56 deletions src/phg/sift/sift.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -109,15 +109,20 @@ std::vector<phg::SIFT::Octave> phg::buildOctaves(const cv::Mat& img, const phg::
// это будет немного быстрее, тк нужно более маленькое ядро свертки на каждый шаг
for (int i = 1; i < n_layers; i++) {
// TODO double sigma_layer = sigma0 * корень из двух нужной степени, чтобы при i==s получали удвоение базового блюра;
double sigma_layer = sigma0 * std::pow(2.0, static_cast<double>(i) / s);
// // вычтем sigma0 чтобы размыть ровно до нужной суммарной сигмы
// TODO sigma_layer = ... (вычитаем как в sigma base);
sigma_layer = std::sqrt(sigma_layer * sigma_layer - sigma0 * sigma0);
// 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);
}

// подготавливаем базовый слой для следующей октавы
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);
base = downsample2x(oct.layers[s]);

// можно использовать и downsample2x_avg(oct.layers[s]), это позволяет потом заапскейлить слои обратно до оригинального разрешения без сдвига
// но потребуется везде изменить формулу для пересчета ключевых точек: pt_upscaled = (pt_downscaled + 0.5) * 2^o - 0.5
Expand All @@ -139,6 +144,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 (size_t 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 @@ -208,16 +216,33 @@ std::vector<cv::KeyPoint> phg::findScaleSpaceExtrema(const std::vector<phg::SIFT
};

// TODO проверить локальный максимум на текущем скейле
for (int i = -1; i <= 1; ++i) {
check(cp[x + i]);
check(cn[x + i]);
if (i != 0) {
check(c[x + i]);
}
}

if (!is_max && !is_min)
continue;

// TODO проверить локальный максимум на предыдущем скейле
for (int i = -1; i <= 1; ++i) {
check(pp[x + i]);
check(pn[x + i]);
check(p[x + i]);
}

if (!is_max && !is_min)
continue;

// TODO проверить локальный максимум на следующем скейле
for (int i = -1; i <= 1; ++i) {
check(np[x + i]);
check(nn[x + i]);
check(n[x + i]);
}

if (!is_max && !is_min)
continue;
Expand All @@ -237,14 +262,21 @@ std::vector<cv::KeyPoint> phg::findScaleSpaceExtrema(const std::vector<phg::SIFT
float ds = (nL.at<float>(yi, xi) - pL.at<float>(yi, xi)) * 0.5f;

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

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 = (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;
float 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 @@ -273,10 +305,10 @@ std::vector<cv::KeyPoint> phg::findScaleSpaceExtrema(const std::vector<phg::SIFT
// из линейной алгебры, сумма диагональных элементов матрицы (след) равна сумме собственных чисел
// определитель матрицы равен произведению собственных чисел
// в случае гессиана (пространственной части: (dxx dxy, dxy, dyy)), собственные числа lambda1, lambda2 - силы кривизны в направлении максимальной кривизны и в ортогональном
// float trace = //TODO ; // = lambda1 + lambda2
// float det = // TODO ; // = lambda1 * lambda2
// if (det <= 0)
// break; // если произведение кривизн отрицательное, то мы находимся в седловой точке, а не в максимуме/минимуме. если нулевое, то это ровная граница вообще
float trace = H(0, 0) + H(1, 1);//TODO ; // = lambda1 + lambda2
float det = H(0, 0) * H(1, 1) - H(0, 1) * H(1, 0);// TODO ; // = lambda1 * lambda2
if (det <= 0)
break; // если произведение кривизн отрицательное, то мы находимся в седловой точке, а не в максимуме/минимуме. если нулевое, то это ровная граница вообще
//
// // если граница незацепистая = грань, то одна кривизна сильно больше чем другая. хотим, чтобы обе кривизны были примерно сопоставимы
// // тогда их отношение r = lambda1/lambda2 будет не очень большим
Expand All @@ -288,6 +320,9 @@ std::vector<cv::KeyPoint> phg::findScaleSpaceExtrema(const std::vector<phg::SIFT
// float r = edge_threshold;
// if (TODO)
// break;
float r = edge_threshold;
if (trace * trace / det > (r + 1) * (r + 1) / r)
break;
}

// скейлим координаты точек обратно до родных размеров картинки
Expand Down Expand Up @@ -379,39 +414,43 @@ 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); // TODO;
float angle = std::atan2(gy, gx); // 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));
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 = TODO;
float bin = angle_deg / 360.0 * 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] += TODO;
// histogram[bin1] += TODO;
histogram[bin0] += (1.0 - frac) * mag * weight;
histogram[bin1] += frac * mag * weight;
}
}

Expand Down Expand Up @@ -450,20 +489,23 @@ std::vector<cv::KeyPoint> phg::computeOrientations(const std::vector<cv::KeyPoin
// 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);
// float offset = TODO;
float a = (left + right - 2 * center) * 0.5;
float b = (right - left) * 0.5;
float offset = -b / (2.0f * 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);
}
}
}
Expand Down Expand Up @@ -579,6 +621,11 @@ std::pair<cv::Mat, std::vector<cv::KeyPoint>> phg::computeDescriptors(const std:
// 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 +656,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 * wy * wx * wo;
}
}
}
Expand All @@ -621,8 +668,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 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