Skip to content

Task02 Леонид Альжанов ITMO #44

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
lalzhanov wants to merge 3 commits into
base: task02
Choose a base branch
from
Open

Task02 Леонид Альжанов ITMO #44

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
7b6eced
tests passing
lalzhanov Mar 2, 2026
9af6a6e
Merge branch 'task01' into task02
lalzhanov Mar 18, 2026
b3b811a
All tests passing
lalzhanov Mar 18, 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
35 changes: 22 additions & 13 deletions src/phg/matching/cl/bruteforce_matcher.cl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,14 +25,18 @@ __kernel void bruteforce_matcher(__global const float* train,
// храним два лучших сопоставления для каждого дескриптора-query:
__local uint res_train_idx_local[KEYPOINTS_PER_WG * 2];
__local float res_distance2_local[KEYPOINTS_PER_WG * 2]; // храним квадраты чтобы не считать корень до самого последнего момента
__local float dist2_for_reduction[NDIM];
// заполняем текущие лучшие дистанции большими значениями
if (dim_id < KEYPOINTS_PER_WG * 2) {
res_distance2_local[dim_id] = FLT_MAX; // полагаемся на то что res_distance2_local размера KEYPOINTS_PER_WG*2==4*2<dim_id<=NDIM==128
}

// грузим 4 дескриптора-query (для каждого из четырех дескрипторов каждый поток грузит значение своей размерности dim_id)
// TODO: т.е. надо прогрузить в query_local все KEYPOINTS_PER_WG=4 дескриптора из query (начиная с индекса query_id0) (а если часть из них выходит за пределы n_query_desc - грузить нули)

for (int query_local_i = 0; query_local_i < KEYPOINTS_PER_WG; query_local_i++) {
query_local[query_local_i * NDIM + dim_id] = query_id0 + query_local_i < n_query_desc ?
query[(query_id0 + query_local_i) * NDIM + dim_id] : 0.0f;
}

barrier(CLK_LOCAL_MEM_FENCE); // дожидаемся прогрузки наших дескрипторов-запросов

for (int train_idx = 0; train_idx < n_train_desc; ++train_idx) {
Expand All @@ -42,14 +46,15 @@ __kernel void bruteforce_matcher(__global const float* train,
// от дескриптора-query в локальной памяти (#query_local_i)
// до дескриптора-train в глобальной памяти (#train_idx)

// TODO посчитать квадрат расстояния по нашей размерности (dim_id) и сохранить его в нашу ячейку в dist2_for_reduction

float d = train_value_dim - query_local[query_local_i * NDIM + dim_id];
dist2_for_reduction[dim_id] = d * d;

barrier(CLK_LOCAL_MEM_FENCE);
// TODO суммируем редукцией все что есть в dist2_for_reduction

int step = NDIM / 2;
while (step > 0) {
if (dim_id < step) {
// TODO
dist2_for_reduction[dim_id] += dist2_for_reduction[dim_id + step];
}
barrier(CLK_LOCAL_MEM_FENCE);
step /= 2;
Expand All @@ -63,17 +68,21 @@ __kernel void bruteforce_matcher(__global const float* train,
#define SECOND_BEST_INDEX 1

// пытаемся улучшить самое лучшее сопоставление для локального дескриптора
if (dist2 <= res_distance2_local[query_local_i * 2 + BEST_INDEX]) {
if (dist2 < res_distance2_local[query_local_i * 2 + BEST_INDEX]) {
// не забываем что прошлое лучшее сопоставление теперь стало вторым по лучшевизне (на данный момент)
res_distance2_local[query_local_i * 2 + SECOND_BEST_INDEX] = res_distance2_local[query_local_i * 2 + BEST_INDEX];
res_train_idx_local[query_local_i * 2 + SECOND_BEST_INDEX] = res_train_idx_local[query_local_i * 2 + BEST_INDEX];
// TODO заменяем нашим (dist2, train_idx) самое лучшее сопоставление для локального дескриптора
} else if (dist2 <= res_distance2_local[query_local_i * 2 + SECOND_BEST_INDEX]) { // может мы улучшили хотя бы второе по лучшевизне сопоставление?
// TODO заменяем второе по лучшевизне сопоставление для локального дескриптора

res_distance2_local[query_local_i * 2 + BEST_INDEX] = dist2;
res_train_idx_local[query_local_i * 2 + BEST_INDEX] = train_idx;
} else if (dist2 < res_distance2_local[query_local_i * 2 + SECOND_BEST_INDEX]) { // может мы улучшили хотя бы второе по лучшевизне сопоставление?
res_distance2_local[query_local_i * 2 + SECOND_BEST_INDEX] = dist2;
res_train_idx_local[query_local_i * 2 + SECOND_BEST_INDEX] = train_idx;
}
}
}
}
barrier(CLK_LOCAL_MEM_FENCE);

// итак, мы нашли два лучших сопоставления для наших KEYPOINTS_PER_WG дескрипторов, надо сохрнить эти результаты в глобальную память
if (dim_id < KEYPOINTS_PER_WG * 2) { // полагаемся на то что нам надо прогрузить KEYPOINTS_PER_WG*2==4*2<dim_id<=NDIM==128
Expand All @@ -82,9 +91,9 @@ __kernel void bruteforce_matcher(__global const float* train,

int query_id = query_id0 + query_local_i;
if (query_id < n_query_desc) {
res_train_idx[query_id * 2 + k] = // TODO
res_query_idx[query_id * 2 + k] = // TODO хм, не масло масленное ли?
res_distance [query_id * 2 + k] = // TODO не забудьте извлечь корень
res_train_idx[query_id * 2 + k] = res_train_idx_local[query_local_i * 2 + k];
res_query_idx[query_id * 2 + k] = query_id;
res_distance [query_id * 2 + k] = sqrt(res_distance2_local[query_local_i * 2 + k]);
}
}
}
109 changes: 70 additions & 39 deletions src/phg/matching/descriptor_matcher.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,15 @@ void phg::DescriptorMatcher::filterMatchesRatioTest(const std::vector<std::vecto
{
filtered_matches.clear();

throw std::runtime_error("not implemented yet");
for (const std::vector<cv::DMatch> &desc_matches : matches) {
if (desc_matches.size() < 2) {
continue;
}

if (desc_matches[0].distance < 0.7 * desc_matches[1].distance) {
filtered_matches.push_back(desc_matches[0]);
}
}
}


Expand All @@ -35,42 +43,65 @@ void phg::DescriptorMatcher::filterMatchesClusters(const std::vector<cv::DMatch>
points_query.at<cv::Point2f>(i) = keypoints_query[matches[i].queryIdx].pt;
points_train.at<cv::Point2f>(i) = keypoints_train[matches[i].trainIdx].pt;
}
//
// // размерность всего 2, так что точное KD-дерево
// std::shared_ptr<cv::flann::IndexParams> index_params = flannKdTreeIndexParams(TODO);
// std::shared_ptr<cv::flann::SearchParams> search_params = flannKsTreeSearchParams(TODO);
//
// std::shared_ptr<cv::flann::Index> index_query = flannKdTreeIndex(points_query, index_params);
// std::shared_ptr<cv::flann::Index> index_train = flannKdTreeIndex(points_train, index_params);
//
// // для каждой точки найти total neighbors ближайших соседей
// cv::Mat indices_query(n_matches, total_neighbours, CV_32SC1);
// cv::Mat distances2_query(n_matches, total_neighbours, CV_32FC1);
// cv::Mat indices_train(n_matches, total_neighbours, CV_32SC1);
// cv::Mat distances2_train(n_matches, total_neighbours, CV_32FC1);
//
// index_query->knnSearch(points_query, indices_query, distances2_query, total_neighbours, *search_params);
// index_train->knnSearch(points_train, indices_train, distances2_train, total_neighbours, *search_params);
//
// // оценить радиус поиска для каждой картинки
// // NB: radius2_query, radius2_train: квадраты радиуса!
// float radius2_query, radius2_train;
// {
// std::vector<double> max_dists2_query(n_matches);
// std::vector<double> max_dists2_train(n_matches);
// for (int i = 0; i < n_matches; ++i) {
// max_dists2_query[i] = distances2_query.at<float>(i, total_neighbours - 1);
// max_dists2_train[i] = distances2_train.at<float>(i, total_neighbours - 1);
// }
//
// int median_pos = n_matches / 2;
// std::nth_element(max_dists2_query.begin(), max_dists2_query.begin() + median_pos, max_dists2_query.end());
// std::nth_element(max_dists2_train.begin(), max_dists2_train.begin() + median_pos, max_dists2_train.end());
//
// radius2_query = max_dists2_query[median_pos] * radius_limit_scale * radius_limit_scale;
// radius2_train = max_dists2_train[median_pos] * radius_limit_scale * radius_limit_scale;
// }
//
// метч остается, если левое и правое множества первых total_neighbors соседей в радиусах поиска(radius2_query, radius2_train) имеют как минимум consistent_matches общих элементов
// // TODO заполнить filtered_matches

// размерность всего 2, так что точное KD-дерево
std::shared_ptr<cv::flann::IndexParams> index_params = flannKdTreeIndexParams(1);
std::shared_ptr<cv::flann::SearchParams> search_params = flannKsTreeSearchParams(std::max(n_matches, 64));

std::shared_ptr<cv::flann::Index> index_query = flannKdTreeIndex(points_query, index_params);
std::shared_ptr<cv::flann::Index> index_train = flannKdTreeIndex(points_train, index_params);

// для каждой точки найти total neighbors ближайших соседей
cv::Mat indices_query(n_matches, total_neighbours, CV_32SC1);
cv::Mat distances2_query(n_matches, total_neighbours, CV_32FC1);
cv::Mat indices_train(n_matches, total_neighbours, CV_32SC1);
cv::Mat distances2_train(n_matches, total_neighbours, CV_32FC1);

index_query->knnSearch(points_query, indices_query, distances2_query, total_neighbours, *search_params);
index_train->knnSearch(points_train, indices_train, distances2_train, total_neighbours, *search_params);

// оценить радиус поиска для каждой картинки
// NB: radius2_query, radius2_train: квадраты радиуса!
float radius2_query, radius2_train;
{
std::vector<double> max_dists2_query(n_matches);
std::vector<double> max_dists2_train(n_matches);
for (int i = 0; i < n_matches; ++i) {
max_dists2_query[i] = distances2_query.at<float>(i, total_neighbours - 1);
max_dists2_train[i] = distances2_train.at<float>(i, total_neighbours - 1);
}

int median_pos = n_matches / 2;
std::nth_element(max_dists2_query.begin(), max_dists2_query.begin() + median_pos, max_dists2_query.end());
std::nth_element(max_dists2_train.begin(), max_dists2_train.begin() + median_pos, max_dists2_train.end());

radius2_query = max_dists2_query[median_pos] * radius_limit_scale * radius_limit_scale;
radius2_train = max_dists2_train[median_pos] * radius_limit_scale * radius_limit_scale;
}

// метч остается, если левое и правое множества первых total_neighbors соседей в радиусах поиска(radius2_query, radius2_train) имеют как минимум consistent_matches общих элементов
for (int i = 0; i < n_matches; i++) {
std::vector <int> query_close;
std::vector <int> train_close;

for (int j = 0; j < total_neighbours; j++) {
if (distances2_query.at<float>(i, j) <= radius2_query) {
query_close.push_back(indices_query.at<int>(i, j));
}
if (distances2_train.at<float>(i, j) <= radius2_train) {
train_close.push_back(indices_train.at<int>(i, j));
}
}

int result = 0;
for (int j = 0; j < train_close.size(); j++) {
if (std::find(query_close.begin(), query_close.end(), train_close[j]) != query_close.end()) {
result++;
}
}

if (result >= consistent_matches) {
filtered_matches.push_back(matches[i]);
}
}
}
18 changes: 15 additions & 3 deletions src/phg/matching/flann_matcher.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,8 @@
phg::FlannMatcher::FlannMatcher()
{
// параметры для приближенного поиска
// index_params = flannKdTreeIndexParams(TODO);
// search_params = flannKsTreeSearchParams(TODO);
index_params = flannKdTreeIndexParams(4);
search_params = flannKsTreeSearchParams(32);
}

void phg::FlannMatcher::train(const cv::Mat &train_desc)
Expand All @@ -17,5 +17,17 @@ void phg::FlannMatcher::train(const cv::Mat &train_desc)

void phg::FlannMatcher::knnMatch(const cv::Mat &query_desc, std::vector<std::vector<cv::DMatch>> &matches, int k) const
{
throw std::runtime_error("not implemented yet");
matches.clear();
int n = query_desc.rows;
cv::Mat inds(n, k, CV_32SC1);
cv::Mat dists(n, k, CV_32FC1);
flann_index -> knnSearch(query_desc, inds, dists, k, *search_params);
matches.resize(n);
for (int i = 0; i < n; i++) {
for (int j = 0; j < k; j++) {
if (int ind = inds.at<int>(i, j); ind > 0) {
matches[i].emplace_back(i, ind, 0, std::sqrt(dists.at<float>(i, j)));
}
}
}
}
114 changes: 60 additions & 54 deletions src/phg/sfm/homography.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,8 +84,8 @@ namespace {
double w1 = ws1[i];

// 8 elements of matrix + free term as needed by gauss routine
// A.push_back({TODO});
// A.push_back({TODO});
A.push_back({0, 0, 0, -x0 * w1, -y0 * w1, -w0 * w1, x0 * y1, y0 * y1, -w0 * y1});
A.push_back({x0 * w1, y0 * w1, w0 * w1, 0, 0, 0, -x0 * x1, -y0 * x1, w0 * x1});
}

int res = gauss(A, H);
Expand Down Expand Up @@ -168,57 +168,57 @@ namespace {
// * (простое описание для понимания)
// * [3] http://ikrisoft.blogspot.com/2015/01/ransac-with-contrario-approach.html

// const int n_matches = points_lhs.size();
//
// // https://en.wikipedia.org/wiki/Random_sample_consensus#Parameters
// const int n_trials = TODO;
//
// const int n_samples = TODO;
// uint64_t seed = 1;
// const double reprojection_error_threshold_px = 2;
//
// int best_support = 0;
// cv::Mat best_H;
//
// std::vector<int> sample;
// for (int i_trial = 0; i_trial < n_trials; ++i_trial) {
// randomSample(sample, n_matches, n_samples, &seed);
//
// cv::Mat H = estimateHomography4Points(points_lhs[sample[0]], points_lhs[sample[1]], points_lhs[sample[2]], points_lhs[sample[3]],
// points_rhs[sample[0]], points_rhs[sample[1]], points_rhs[sample[2]], points_rhs[sample[3]]);
//
// int support = 0;
// for (int i_point = 0; i_point < n_matches; ++i_point) {
// try {
// cv::Point2d proj = phg::transformPoint(points_lhs[i_point], H);
// if (cv::norm(proj - cv::Point2d(points_rhs[i_point])) < reprojection_error_threshold_px) {
// ++support;
// }
// } catch (const std::exception &e)
// {
// std::cerr << e.what() << std::endl;
// }
// }
//
// if (support > best_support) {
// best_support = support;
// best_H = H;
//
// std::cout << "estimateHomographyRANSAC : support: " << best_support << "/" << n_matches << std::endl;
//
// if (best_support == n_matches) {
// break;
// }
// }
// }
//
// std::cout << "estimateHomographyRANSAC : best support: " << best_support << "/" << n_matches << std::endl;
//
// if (best_support == 0) {
// throw std::runtime_error("estimateHomographyRANSAC : failed to estimate homography");
// }
//
// return best_H;
const int n_matches = points_lhs.size();

// https://en.wikipedia.org/wiki/Random_sample_consensus#Parameters
const int n_trials = 1000;

const int n_samples = 5;
uint64_t seed = 1;
const double reprojection_error_threshold_px = 2;

int best_support = 0;
cv::Mat best_H;

std::vector<int> sample;
for (int i_trial = 0; i_trial < n_trials; ++i_trial) {
randomSample(sample, n_matches, n_samples, &seed);

cv::Mat H = estimateHomography4Points(points_lhs[sample[0]], points_lhs[sample[1]], points_lhs[sample[2]], points_lhs[sample[3]],
points_rhs[sample[0]], points_rhs[sample[1]], points_rhs[sample[2]], points_rhs[sample[3]]);

int support = 0;
for (int i_point = 0; i_point < n_matches; ++i_point) {
try {
cv::Point2d proj = phg::transformPoint(points_lhs[i_point], H);
if (cv::norm(proj - cv::Point2d(points_rhs[i_point])) < reprojection_error_threshold_px) {
++support;
}
} catch (const std::exception &e)
{
std::cerr << e.what() << std::endl;
}
}

if (support > best_support) {
best_support = support;
best_H = H;

std::cout << "estimateHomographyRANSAC : support: " << best_support << "/" << n_matches << std::endl;

if (best_support == n_matches) {
break;
}
}
}

std::cout << "estimateHomographyRANSAC : best support: " << best_support << "/" << n_matches << std::endl;

if (best_support == 0) {
throw std::runtime_error("estimateHomographyRANSAC : failed to estimate homography");
}

return best_H;
}

}
Expand All @@ -238,7 +238,13 @@ cv::Mat phg::findHomographyCV(const std::vector<cv::Point2f> &points_lhs, const
// таким преобразованием внутри занимается функции cv::perspectiveTransform и cv::warpPerspective
cv::Point2d phg::transformPoint(const cv::Point2d &pt, const cv::Mat &T)
{
throw std::runtime_error("not implemented yet");
const double x = T.at<double>(0, 0) * pt.x + T.at<double>(0, 1) * pt.y + T.at<double>(0, 2);
const double y = T.at<double>(1, 0) * pt.x + T.at<double>(1, 1) * pt.y + T.at<double>(1, 2);
const double z = T.at<double>(2, 0) * pt.x + T.at<double>(2, 1) * pt.y + T.at<double>(2, 2);
if (std::abs(z) < 1e-12) {
throw std::runtime_error("phg::transformPoint: scale too small");
}
return {x / z, y / z};
}

cv::Point2d phg::transformPointCV(const cv::Point2d &pt, const cv::Mat &T) {
Expand Down
Loading
Loading