Skip to content

Task01 Петр Маяцкий ИТМО #5

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
BellatorHonoris07 wants to merge 1 commit into
base: task01
Choose a base branch
from
Open

Task01 Петр Маяцкий ИТМО #5

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
197 changes: 97 additions & 100 deletions src/phg/sift/sift.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,21 +107,17 @@ std::vector<phg::SIFT::Octave> phg::buildOctaves(const cv::Mat& img, const phg::
// для простоты в каждой октаве будем каждый раз блюрить базовую картинку с полной сигмой
// можно подумать, как сделать эффективнее - для построения n+1 слоя доблюревать уже поблюренный n-ый слой, так чтобы в итоге получилась такая же сигма
// это будет немного быстрее, тк нужно более маленькое ядро свертки на каждый шаг
double k = std::pow(2, (double) 1 / s);
double sigma_prev = sigma0;
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 = std::sqrt(sigma_prev * sigma_prev * k * k - sigma_prev * sigma_prev);
cv::GaussianBlur(oct.layers[i-1], oct.layers[i], cv::Size(), sigma_layer, sigma_layer);
sigma_prev *= k;
}

// подготавливаем базовый слой для следующей октавы
if (o + 1 < n_octaves) {
// используется в opencv, формула для пересчета ключевых точек: pt_upscaled = 2^o * pt_downscaled
// TODO cv::resize(даунскейлим текущий слой в два раза, без интерполяции, просто сабсепмлинг);

// можно использовать и downsample2x_avg(oct.layers[s]), это позволяет потом заапскейлить слои обратно до оригинального разрешения без сдвига
// но потребуется везде изменить формулу для пересчета ключевых точек: pt_upscaled = (pt_downscaled + 0.5) * 2^o - 0.5

cv::resize(oct.layers[s], base, cv::Size(oct.layers[s].cols/2, oct.layers[s].rows/2), 0, 0, cv::INTER_NEAREST);
if (verbose_level)
std::cout << "new octave base size: " << base.size().width << std::endl;
}
Expand All @@ -133,12 +129,12 @@ std::vector<phg::SIFT::Octave> phg::buildOctaves(const cv::Mat& img, const phg::
std::vector<phg::SIFT::Octave> phg::buildDoG(const std::vector<phg::SIFT::Octave>& octaves, const phg::SIFTParams& p, int verbose_level)
{
std::vector<phg::SIFT::Octave> 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);

// 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 @@ -207,17 +203,28 @@ std::vector<cv::KeyPoint> phg::findScaleSpaceExtrema(const std::vector<phg::SIFT
is_min = false;
};

// TODO проверить локальный максимум на текущем скейле
auto check_rows = [&](const float* row_up, const float* row_mid, const float* row_down, bool skip_center) {
for (int i = -1; i <= 1; i++) {
check(row_up[x + i]);
check(row_down[x + i]);
if (skip_center && i == 0) {
continue;
}
check(row_mid[x + i]);
}
};

check_rows(cp, c, cn, true);

if (!is_max && !is_min)
continue;

// TODO проверить локальный максимум на предыдущем скейле
check_rows(pp, p, pn, false);

if (!is_max && !is_min)
continue;

// TODO проверить локальный максимум на следующем скейле
check_rows(np, n, nn, false);

if (!is_max && !is_min)
continue;
Expand All @@ -237,14 +244,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 All @@ -270,26 +276,18 @@ std::vector<cv::KeyPoint> phg::findScaleSpaceExtrema(const std::vector<phg::SIFT

// фильтрация по зацепистости
if (params.enable_edge_like_filtering) {
// из линейной алгебры, сумма диагональных элементов матрицы (след) равна сумме собственных чисел
// определитель матрицы равен произведению собственных чисел
// в случае гессиана (пространственной части: (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;
float det = dxx * dyy - dxy * dxy;

if (det <= 0) {
break;
}

float r = edge_threshold;
if ((trace*trace)/det >= ((r + 1) * (r + 1) / r)) {
break;
}
}
// скейлим координаты точек обратно до родных размеров картинки
// !!! если выбираем при даунскейле другую политику, с усреднением вместо ресемплинга, то надо здесь применять формулу со сдвигами на полпикселя
float scale = (real_octave >= 0) ? (float)(1 << real_octave) : (1.f / (float)(1 << (-real_octave)));
Expand Down Expand Up @@ -379,39 +377,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;
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 +448,20 @@ 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 = (left-right)/(2*(left-2*center+right));
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 +572,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 +607,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 @@ -620,9 +618,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
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