Skip to content

mvp #8

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
VashchilkoAV wants to merge 1 commit into
base: task01
Choose a base branch
from
+108 −80
Open

mvp #8

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
186 changes: 107 additions & 79 deletions src/phg/sift/sift.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -108,14 +108,18 @@ 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 получали удвоение базового блюра;
double sigma_layer = sigma0 * std::pow((double) 2, (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);
// double sigma_layer = sigma0 * корень из двух нужной степени, чтобы при i==s получали удвоение базового блюра;
// // вычтем sigma0 чтобы размыть ровно до нужной суммарной сигмы
// TODO sigma_layer = ... (вычитаем как в sigma base);
// sigma_layer = ... (вычитаем как в sigma base);
// cv::GaussianBlur(oct.layers[0], oct.layers[i], cv::Size(), sigma_layer, sigma_layer);
}

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

Expand All @@ -137,8 +141,10 @@ 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);

// 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];
}
// каждый слой дога это разница n+1 и n-й гауссианы
}

return dog;
Expand Down Expand Up @@ -207,17 +213,40 @@ std::vector<cv::KeyPoint> phg::findScaleSpaceExtrema(const std::vector<phg::SIFT
is_min = false;
};

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

check(c[x - 1]);
check(c[x + 1]);
check(cn[x - 1]);
check(cn[x]);
check(cn[x + 1]);
check(cp[x - 1]);
check(cp[x]);
check(cp[x + 1]);

if (!is_max && !is_min)
continue;

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

if (!is_max && !is_min)
continue;

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

if (!is_max && !is_min)
continue;
Expand All @@ -238,13 +267,13 @@ 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 = pL.at<float>(yi, xi) + nL.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 @@ -273,10 +302,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 = dxx + dyy; // = lambda1 + lambda2
float det = dxx * dyy - dxy * dxy; // = lambda1 * lambda2
if (det <= 0)
break; // если произведение кривизн отрицательное, то мы находимся в седловой точке, а не в максимуме/минимуме. если нулевое, то это ровная граница вообще
//
// // если граница незацепистая = грань, то одна кривизна сильно больше чем другая. хотим, чтобы обе кривизны были примерно сопоставимы
// // тогда их отношение r = lambda1/lambda2 будет не очень большим
Expand All @@ -285,9 +314,9 @@ std::vector<cv::KeyPoint> phg::findScaleSpaceExtrema(const std::vector<phg::SIFT
// // и просто как интуиция, при больших r это выражение просто до r сокращается
//
// // в итоге получается что порог edge_threshold в отличие от response_threshold наоборот, чем больше тем расслабленнее
// float r = edge_threshold;
// if (TODO)
// break;
float r = edge_threshold;
if (trace * trace / det > ((double) r + 1) * ((double) r + 1) / r)
break;
}

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

Expand Down Expand Up @@ -450,20 +479,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) / (left + right - 2.f * center) / 2.f;
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 +603,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 +638,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] += wy * wx * wo * weighted_mag;
}
}
}
Expand All @@ -620,9 +649,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