Skip to content

Task01 Леонид Тернопол СПбГУ #6

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
firelion9 wants to merge 1 commit into from
Closed

Task01 Леонид Тернопол СПбГУ #6

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
3 changes: 3 additions & 0 deletions .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,3 +2,6 @@
.idea
build
cmake-build*
/.cache/
/.venv/
/.vscode
211 changes: 112 additions & 99 deletions src/phg/sift/sift.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#include "sift.h"
#include "libutils/rasserts.h"

#include <cmath>
#include <iostream>
#include <numeric>
#include <opencv2/features2d.hpp>
Expand Down Expand Up @@ -108,16 +109,16 @@ std::vector<phg::SIFT::Octave> 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 sigma_layer = sigma0 * std::pow(2.0, (double) i / s);
// вычтем sigma0 чтобы размыть ровно до нужной суммарной сигмы
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, oct.layers[0].size() / 2, 0, 0, cv::INTER_NEAREST);

// можно использовать и downsample2x_avg(oct.layers[s]), это позволяет потом заапскейлить слои обратно до оригинального разрешения без сдвига
// но потребуется везде изменить формулу для пересчета ключевых точек: pt_upscaled = (pt_downscaled + 0.5) * 2^o - 0.5
Expand All @@ -136,9 +137,13 @@ std::vector<phg::SIFT::Octave> phg::buildDoG(const std::vector<phg::SIFT::Octave

for (size_t o = 0; o < octaves.size(); o++) {
const phg::SIFT::Octave& octave = octaves[o];
dog[o].layers.resize(octave.layers.size() - 1);
phg::SIFT::Octave& dg = dog[o];
dg.layers.resize(octave.layers.size() - 1);

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

return dog;
Expand Down Expand Up @@ -206,21 +211,28 @@ std::vector<cv::KeyPoint> phg::findScaleSpaceExtrema(const std::vector<phg::SIFT
if (v <= val)
is_min = false;
};
auto check_layer = [&](const float *y1, const float *y0, const float *y2) {
is_max = true;
is_min = true;

check(y0[x - 1]);
check(y0[x]);
check(y0[x + 1]);
check(y1[x - 1]);
check(y1[x + 1]);
check(y2[x - 1]);
check(y2[x]);
check(y2[x + 1]);
};

// TODO проверить локальный максимум на текущем скейле

if (!is_max && !is_min)
continue;

// TODO проверить локальный максимум на предыдущем скейле

if (!is_max && !is_min)
continue;
check_layer(c, cp, cn);
if (!is_max && !is_min) continue;

// TODO проверить локальный максимум на следующем скейле
check_layer(p, pp, pn);
if (!is_max && !is_min) continue;

if (!is_max && !is_min)
continue;
check_layer(n, np, nn);
if (!is_max && !is_min) continue;

int xi = x, yi = y, li = layer;

Expand All @@ -237,14 +249,13 @@ 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 = 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,21 +284,22 @@ 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; // если произведение кривизн отрицательное, то мы находимся в седловой точке, а не в максимуме/минимуме. если нулевое, то это ровная граница вообще
//
// // если граница незацепистая = грань, то одна кривизна сильно больше чем другая. хотим, чтобы обе кривизны были примерно сопоставимы
// // тогда их отношение r = lambda1/lambda2 будет не очень большим
// // если расписать trace * trace / det через r, то получится (r + 1) ^ 2 / r
// // функция растущая по r, так что если наше фактическое значение trace * trace / det выше (r + 1) ^ 2 / r, то и наше отношение кривизн больше порога, значит плохая зацепистость
// // и просто как интуиция, при больших r это выражение просто до r сокращается
//
// // в итоге получается что порог edge_threshold в отличие от response_threshold наоборот, чем больше тем расслабленнее
// float r = edge_threshold;
// if (TODO)
// break;
float trace = dxx + dyy; // = lambda1 + lambda2
float det = dxx * dyy - dxy * dxy; // = lambda1 * lambda2
if (det <= 0)
break; // если произведение кривизн отрицательное, то мы находимся в седловой точке, а не в максимуме/минимуме. если нулевое, то это ровная граница вообще

// если граница незацепистая = грань, то одна кривизна сильно больше чем другая. хотим, чтобы обе кривизны были примерно сопоставимы
// тогда их отношение r = lambda1/lambda2 будет не очень большим
// если расписать trace * trace / det через r, то получится (r + 1) ^ 2 / r
// функция растущая по r, так что если наше фактическое значение trace * trace / det выше (r + 1) ^ 2 / r, то и наше отношение кривизн больше порога, значит плохая зацепистость
// и просто как интуиция, при больших r это выражение просто до r сокращается

// в итоге получается что порог edge_threshold в отличие от response_threshold наоборот, чем больше тем расслабленнее
float r = trace * trace / det;
if (r < 1.0) r = 1 / r;
if (r > edge_threshold)
break;
}

// скейлим координаты точек обратно до родных размеров картинки
Expand Down Expand Up @@ -379,39 +391,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::hypot(gx, 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 / 360.0f * 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] += (1 - frac) * weight * mag;
histogram[bin1] += frac * weight * mag;
}
}

Expand Down Expand Up @@ -450,20 +462,22 @@ 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 a = (left + right - 2 * center) / 2;
float b = (right - left) / 2;
float offset = -b / (2 * 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 @@ -574,11 +588,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 +623,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 @@ -620,9 +634,8 @@ std::pair<cv::Mat, std::vector<cv::KeyPoint>> 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;
}
}
}
Expand Down
5 changes: 1 addition & 4 deletions tests/test_sift.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,10 +25,7 @@

#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

Expand Down
Loading