Skip to content

Commit 4e4e451

Browse files
HackAss2002HackAss2002
committed
Task01 Вячеслав Григорович ITMO
1 parent 91387d1 commit 4e4e451

2 files changed

Lines changed: 122 additions & 92 deletions

File tree

‎src/phg/sift/sift.cpp‎

Lines changed: 120 additions & 88 deletions
Original file line numberDiff line numberDiff line change
@@ -108,16 +108,16 @@ std::vector<phg::SIFT::Octave> phg::buildOctaves(const cv::Mat& img, const phg::
108108
// можно подумать, как сделать эффективнее - для построения n+1 слоя доблюревать уже поблюренный n-ый слой, так чтобы в итоге получилась такая же сигма
109109
// это будет немного быстрее, тк нужно более маленькое ядро свертки на каждый шаг
110110
for (int i = 1; i < n_layers; i++) {
111-
// TODO double sigma_layer = sigma0 * корень из двух нужной степени, чтобы при i==s получали удвоение базового блюра;
111+
double sigma_layer = sigma0 * std::pow(2.0, static_cast<double>(i) / s);
112112
// // вычтем sigma0 чтобы размыть ровно до нужной суммарной сигмы
113-
// TODO sigma_layer = ... (вычитаем как в sigma base);
114-
// cv::GaussianBlur(oct.layers[0], oct.layers[i], cv::Size(), sigma_layer, sigma_layer);
113+
sigma_layer = std::sqrt(sigma_layer * sigma_layer - sigma0 * sigma0);
114+
cv::GaussianBlur(oct.layers[0], oct.layers[i], cv::Size(), sigma_layer, sigma_layer);
115115
}
116116

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

122122
// можно использовать и downsample2x_avg(oct.layers[s]), это позволяет потом заапскейлить слои обратно до оригинального разрешения без сдвига
123123
// но потребуется везде изменить формулу для пересчета ключевых точек: pt_upscaled = (pt_downscaled + 0.5) * 2^o - 0.5
@@ -138,7 +138,9 @@ std::vector<phg::SIFT::Octave> phg::buildDoG(const std::vector<phg::SIFT::Octave
138138
const phg::SIFT::Octave& octave = octaves[o];
139139
dog[o].layers.resize(octave.layers.size() - 1);
140140

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

144146
return dog;
@@ -207,17 +209,49 @@ std::vector<cv::KeyPoint> phg::findScaleSpaceExtrema(const std::vector<phg::SIFT
207209
is_min = false;
208210
};
209211

210-
// TODO проверить локальный максимум на текущем скейле
212+
// проверить локальный максимум на текущем скейле
213+
check(cp[x - 1]);
214+
check(cp[x]);
215+
check(cp[x + 1]);
216+
217+
check(c[x - 1]);
218+
check(c[x + 1]);
219+
220+
check(cn[x - 1]);
221+
check(cn[x]);
222+
check(cn[x + 1]);
211223

212224
if (!is_max && !is_min)
213225
continue;
214226

215-
// TODO проверить локальный максимум на предыдущем скейле
227+
// проверить локальный максимум на предыдущем скейле
228+
check(pp[x - 1]);
229+
check(pp[x]);
230+
check(pp[x + 1]);
231+
232+
check(p[x - 1]);
233+
check(p[x]);
234+
check(p[x + 1]);
235+
236+
check(pn[x - 1]);
237+
check(pn[x]);
238+
check(pn[x + 1]);
216239

217240
if (!is_max && !is_min)
218241
continue;
219242

220-
// TODO проверить локальный максимум на следующем скейле
243+
// проверить локальный максимум на следующем скейле
244+
check(np[x - 1]);
245+
check(np[x]);
246+
check(np[x + 1]);
247+
248+
check(n[x - 1]);
249+
check(n[x]);
250+
check(n[x + 1]);
251+
252+
check(nn[x - 1]);
253+
check(nn[x]);
254+
check(nn[x + 1]);
221255

222256
if (!is_max && !is_min)
223257
continue;
@@ -237,14 +271,13 @@ std::vector<cv::KeyPoint> phg::findScaleSpaceExtrema(const std::vector<phg::SIFT
237271
float ds = (nL.at<float>(yi, xi) - pL.at<float>(yi, xi)) * 0.5f;
238272

239273
// гессиан
240-
float dxx, dxy, dyy, dxs, dys, dss;
241-
// float dxx = cL.at<float>(yi, xi + 1) + cL.at<float>(yi, xi - 1) - 2.f * resp_center;
242-
// float dyy = TODO;
243-
// float dss = TODO;
244-
//
245-
// 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;
246-
// float dxs = TODO;
247-
// float dys = TODO;
274+
float dxx = cL.at<float>(yi, xi + 1) + cL.at<float>(yi, xi - 1) - 2.f * resp_center;
275+
float dyy = cL.at<float>(yi + 1, xi) + cL.at<float>(yi - 1, xi) - 2.f * resp_center;
276+
float dss = nL.at<float>(yi, xi) + pL.at<float>(yi, xi) - 2.f * resp_center;
277+
278+
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;
279+
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;
280+
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;
248281

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

@@ -273,21 +306,21 @@ std::vector<cv::KeyPoint> phg::findScaleSpaceExtrema(const std::vector<phg::SIFT
273306
// из линейной алгебры, сумма диагональных элементов матрицы (след) равна сумме собственных чисел
274307
// определитель матрицы равен произведению собственных чисел
275308
// в случае гессиана (пространственной части: (dxx dxy, dxy, dyy)), собственные числа lambda1, lambda2 - силы кривизны в направлении максимальной кривизны и в ортогональном
276-
// float trace = //TODO ; // = lambda1 + lambda2
277-
// float det = // TODO ; // = lambda1 * lambda2
278-
// if (det <= 0)
279-
// break; // если произведение кривизн отрицательное, то мы находимся в седловой точке, а не в максимуме/минимуме. если нулевое, то это ровная граница вообще
280-
//
281-
// // если граница незацепистая = грань, то одна кривизна сильно больше чем другая. хотим, чтобы обе кривизны были примерно сопоставимы
282-
// // тогда их отношение r = lambda1/lambda2 будет не очень большим
283-
// // если расписать trace * trace / det через r, то получится (r + 1) ^ 2 / r
284-
// // функция растущая по r, так что если наше фактическое значение trace * trace / det выше (r + 1) ^ 2 / r, то и наше отношение кривизн больше порога, значит плохая зацепистость
285-
// // и просто как интуиция, при больших r это выражение просто до r сокращается
286-
//
287-
// // в итоге получается что порог edge_threshold в отличие от response_threshold наоборот, чем больше тем расслабленнее
288-
// float r = edge_threshold;
289-
// if (TODO)
290-
// break;
309+
float trace = dxx + dyy ; // = lambda1 + lambda2
310+
float det = dxx * dyy - dxy * dxy ; // = lambda1 * lambda2
311+
if (det <= 0)
312+
break; // если произведение кривизн отрицательное, то мы находимся в седловой точке, а не в максимуме/минимуме. если нулевое, то это ровная граница вообще
313+
314+
// если граница незацепистая = грань, то одна кривизна сильно больше чем другая. хотим, чтобы обе кривизны были примерно сопоставимы
315+
// тогда их отношение r = lambda1/lambda2 будет не очень большим
316+
// если расписать trace * trace / det через r, то получится (r + 1) ^ 2 / r
317+
// функция растущая по r, так что если наше фактическое значение trace * trace / det выше (r + 1) ^ 2 / r, то и наше отношение кривизн больше порога, значит плохая зацепистость
318+
// и просто как интуиция, при больших r это выражение просто до r сокращается
319+
320+
// в итоге получается что порог edge_threshold в отличие от response_threshold наоборот, чем больше тем расслабленнее
321+
float r = edge_threshold;
322+
if (trace * trace / det > (r + 1) * (r + 1) / r)
323+
break;
291324
}
292325

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

380413
for (int dy = -radius; dy <= radius; dy++) {
381414
for (int dx = -radius; dx <= radius; dx++) {
382-
// int px = xi + dx;
383-
// int py = yi + dy;
384-
//
385-
// // градиент
386-
// float gx = img.at<float>(py, px + 1) - img.at<float>(py, px - 1);
387-
// float gy = img.at<float>(py + 1, px) - img.at<float>(py - 1, px);
388-
//
389-
// float mag = TODO;
390-
// float angle = std::atan2(TODO); // [-pi, pi]
391-
//
392-
// float angle_deg = angle * 180.f / (float) CV_PI;
393-
// if (angle_deg < 0.f) angle_deg += 360.f;
394-
//
395-
// // гауссово взвешивание голоса точки с затуханием к краям
396-
// float weight = std::exp(-(TODO) / (2.f * sigma_win * sigma_win));
397-
// if (!params.enable_orientation_gaussian_weighting) {
398-
// weight = 1.f;
399-
// }
400-
//
401-
// // голосуем в гистограмме направлений. находим два ближайших бина и гладко распределяем голос между ними
402-
// // в таком случае, голос попавший близко к границе между бинами, проголосует поровну за оба бина
403-
// float bin = TODO;
404-
// if (bin >= n_bins) bin -= n_bins;
405-
// int bin0 = (int) bin;
406-
// int bin1 = (bin0 + 1) % n_bins;
407-
//
408-
// float frac = bin - bin0;
409-
// if (!params.enable_orientation_bin_interpolation) {
410-
// frac = 0.f;
411-
// }
412-
//
413-
// histogram[bin0] += TODO;
414-
// histogram[bin1] += TODO;
415+
int px = xi + dx;
416+
int py = yi + dy;
417+
418+
// градиент
419+
float gx = img.at<float>(py, px + 1) - img.at<float>(py, px - 1);
420+
float gy = img.at<float>(py + 1, px) - img.at<float>(py - 1, px);
421+
422+
float mag = std::sqrt(gx * gx + gy * gy);
423+
float angle = std::atan2(gy, gx); // [-pi, pi]
424+
425+
float angle_deg = angle * 180.f / (float) CV_PI;
426+
if (angle_deg < 0.f) angle_deg += 360.f;
427+
428+
// гауссово взвешивание голоса точки с затуханием к краям
429+
float weight = std::exp(-(dx * dx + dy * dy) / (2.f * sigma_win * sigma_win));
430+
if (!params.enable_orientation_gaussian_weighting) {
431+
weight = 1.f;
432+
}
433+
434+
// голосуем в гистограмме направлений. находим два ближайших бина и гладко распределяем голос между ними
435+
// в таком случае, голос попавший близко к границе между бинами, проголосует поровну за оба бина
436+
float bin = n_bins * angle_deg / 360.0;
437+
if (bin >= n_bins) bin -= n_bins;
438+
int bin0 = (int) bin;
439+
int bin1 = (bin0 + 1) % n_bins;
440+
441+
float frac = bin - bin0;
442+
if (!params.enable_orientation_bin_interpolation) {
443+
frac = 0.f;
444+
}
445+
446+
histogram[bin0] += (1.0 - frac) * mag * weight;
447+
histogram[bin1] += frac * mag * weight;
415448
}
416449
}
417450

@@ -450,20 +483,20 @@ std::vector<cv::KeyPoint> phg::computeOrientations(const std::vector<cv::KeyPoin
450483
// f(1) + f(-1) = 2a + 2c -> a = (left + right - 2 * center) / 2
451484
// f(1) - f(-1) = 2b -> b = (right - left) / 2
452485

453-
// float offset = TODO;
454-
// if (!params.enable_orientation_subpixel_localization) {
455-
// offset = 0.f;
456-
// }
457-
//
458-
// float bin_real = i + offset;
459-
// if (bin_real < 0.f) bin_real += n_bins;
460-
// if (bin_real >= n_bins) bin_real -= n_bins;
461-
//
462-
// float angle = bin_real * 360.f / n_bins;
463-
//
464-
// cv::KeyPoint new_kp = kp;
465-
// new_kp.angle = angle;
466-
// oriented_kpts.push_back(new_kp);
486+
float offset = -(right - left) / (2 * (left + right - 2 * center));
487+
if (!params.enable_orientation_subpixel_localization) {
488+
offset = 0.f;
489+
}
490+
491+
float bin_real = i + offset;
492+
if (bin_real < 0.f) bin_real += n_bins;
493+
if (bin_real >= n_bins) bin_real -= n_bins;
494+
495+
float angle = bin_real * 360.f / n_bins;
496+
497+
cv::KeyPoint new_kp = kp;
498+
new_kp.angle = angle;
499+
oriented_kpts.push_back(new_kp);
467500
}
468501
}
469502
}
@@ -574,11 +607,11 @@ std::pair<cv::Mat, std::vector<cv::KeyPoint>> phg::computeDescriptors(const std:
574607
bin_o -= n_orient_bins;
575608

576609
// семплы вблизи края патча взвешиваем с меньшим весом
577-
// float weight = std::exp(-(TODO) / (2.f * sigma_desc * sigma_desc));
578-
// if (!params.enable_descriptor_gaussian_weighting) {
579-
// weight = 1.f;
580-
// }
581-
// float weighted_mag = mag * weight;
610+
float weight = std::exp(-(rot_x * rot_x + rot_y * rot_y) / (2.f * sigma_desc * sigma_desc));
611+
if (!params.enable_descriptor_gaussian_weighting) {
612+
weight = 1.f;
613+
}
614+
float weighted_mag = mag * weight;
582615

583616
if (params.enable_descriptor_bin_interpolation) {
584617
// размажем вклад weighted_mag по пространственным бинам и по бинам гистограммок трилинейной интерполяцией
@@ -609,8 +642,8 @@ std::pair<cv::Mat, std::vector<cv::KeyPoint>> phg::computeDescriptors(const std:
609642
io += n_orient_bins;
610643
float wo = (dio == 0) ? (1.f - fo) : fo;
611644

612-
// int idx = TODO;
613-
// desc[idx] += TODO;
645+
int idx = (iy * n_spatial_bins + ix) * n_orient_bins + io;
646+
desc[idx] += weighted_mag * wx * wy * wo;
614647
}
615648
}
616649
}
@@ -620,9 +653,8 @@ std::pair<cv::Mat, std::vector<cv::KeyPoint>> phg::computeDescriptors(const std:
620653
int io_nearest = (int)std::round(bin_o) % n_orient_bins;
621654

622655
if (ix_nearest >= 0 && ix_nearest < n_spatial_bins && iy_nearest >= 0 && iy_nearest < n_spatial_bins) {
623-
// TODO uncomment
624-
// int idx = (iy_nearest * n_spatial_bins + ix_nearest) * n_orient_bins + io_nearest;
625-
// desc[idx] += weighted_mag;
656+
int idx = (iy_nearest * n_spatial_bins + ix_nearest) * n_orient_bins + io_nearest;
657+
desc[idx] += weighted_mag;
626658
}
627659
}
628660
}

‎tests/test_sift.cpp‎

Lines changed: 2 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -25,10 +25,7 @@
2525

2626
#define GAUSSIAN_NOISE_STDDEV 1.0
2727

28-
// TODO ENABLE ME
29-
// TODO ENABLE ME
30-
// TODO ENABLE ME
31-
#define ENABLE_MY_SIFT_TESTING 0
28+
#define ENABLE_MY_SIFT_TESTING 1
3229

3330
#define DENY_CREATE_REF_DATA 1
3431

@@ -914,6 +911,7 @@ TEST(SIFT, PairMatching)
914911

915912
phg::SIFTParams params;
916913
params.nfeatures = 10000;
914+
params.orient_peak_ratio = 0.625;
917915

918916
std::cout << "matching using opencv orb..." << std::endl;
919917
auto orb_cv = cv::ORB::create(params.nfeatures);

0 commit comments

Comments
 (0)