Skip to content

Task02 Илья Власов ИТМО #40

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

Task02 Илья Власов ИТМО #40

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
aee8083
+hw1
Malchicisinternetov Mar 1, 2026
ceef316
Merge branch 'task01' into task02
Malchicisinternetov Mar 14, 2026
ab4a3fd
+hw2
Malchicisinternetov Mar 17, 2026
bc561b4
fix gpu
Malchicisinternetov Mar 17, 2026
8573f23
fix gpu 2
Malchicisinternetov Mar 17, 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
34 changes: 22 additions & 12 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 @@ -28,11 +28,19 @@ __kernel void bruteforce_matcher(__global const float* train,
// заполняем текущие лучшие дистанции большими значениями
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
res_train_idx_local[dim_id] = 0;
}

__local float dist2_for_reduction[NDIM];

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

for(int i = 0; i < KEYPOINTS_PER_WG; i++) {
if (query_id0 + i < n_query_desc)
query_local[i * NDIM + dim_id] = query[(query_id0 + i) * NDIM + dim_id];
else
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 +50,14 @@ __kernel void bruteforce_matcher(__global const float* train,
// от дескриптора-query в локальной памяти (#query_local_i)
// до дескриптора-train в глобальной памяти (#train_idx)

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

float distance2 = train_value_dim - query_local[query_local_i * NDIM + dim_id];
dist2_for_reduction[dim_id] = distance2 * distance2;
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,13 +71,15 @@ __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;
}
}
}
Expand All @@ -82,9 +92,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]);
}
}
}
115 changes: 74 additions & 41 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,12 @@ void phg::DescriptorMatcher::filterMatchesRatioTest(const std::vector<std::vecto
{
filtered_matches.clear();

throw std::runtime_error("not implemented yet");
for (const auto &match : matches) {
if (match.size() < 2) continue;
if (constexpr double ratio = 0.7; match[0].distance < ratio*match[1].distance) {
filtered_matches.push_back(match[0]);
}
}
}


Expand All @@ -32,45 +37,73 @@ void phg::DescriptorMatcher::filterMatchesClusters(const std::vector<cv::DMatch>
cv::Mat points_query(n_matches, 2, CV_32FC1);
cv::Mat points_train(n_matches, 2, CV_32FC1);
for (int i = 0; i < n_matches; ++i) {
points_query.at<cv::Point2f>(i) = keypoints_query[matches[i].queryIdx].pt;
points_train.at<cv::Point2f>(i) = keypoints_train[matches[i].trainIdx].pt;
points_query.at<float>(i, 0) = keypoints_query[matches[i].queryIdx].pt.x;
points_query.at<float>(i, 1) = keypoints_query[matches[i].queryIdx].pt.y;
points_train.at<float>(i, 0) = keypoints_train[matches[i].trainIdx].pt.x;
points_train.at<float>(i, 1) = keypoints_train[matches[i].trainIdx].pt.y;
}

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

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));
}
}

for (int j = 0; j < total_neighbours; 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]);
}
}
//
// // размерность всего 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
}
23 changes: 20 additions & 3 deletions src/phg/matching/flann_matcher.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,13 +1,14 @@
#include <iostream>
#include "flann_matcher.h"
#include "flann_factory.h"
#include "gms_matcher_impl.h"


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 +18,21 @@ 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 indices_query(n, k, CV_32SC1);
cv::Mat distances_query(n, k, CV_32FC1);
flann_index->knnSearch(query_desc, indices_query, distances_query, k, *search_params);
matches.resize(n);
for (int i = 0; i < n; i++) {
matches.reserve(k);
for (int j = 0; j < k; j++) {
int ind = indices_query.at<int>(i, j);
if (ind < 0) continue;
float dist = std::sqrt(distances_query.at<float>(i, j));
matches[i].emplace_back(i, ind, 0, dist);

}
//std::sort(matches[i].begin(), matches[i].end(), [](const cv::DMatch a, const cv::DMatch b) {return a.distance < b.distance;});
}
}
120 changes: 65 additions & 55 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 @@ -162,63 +162,65 @@ namespace {
throw std::runtime_error("findHomography: points_lhs.size() != points_rhs.size()");
}

// TODO Дополнительный балл, если вместо обычной версии будет использована модификация a-contrario RANSAC
// * [1] Automatic Homographic Registration of a Pair of Images, with A Contrario Elimination of Outliers. (Lionel Moisan, Pierre Moulon, Pascal Monasse)
// * [2] Adaptive Structure from Motion with a contrario model estimation. (Pierre Moulon, Pascal Monasse, Renaud Marlet)
// * (простое описание для понимания)
// * [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 = 4;
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) {
try {
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;
}
}
} catch (...) {
}


}

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 +240,15 @@ 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");
std::vector <double> res(3);
res[0] = T.at<double>(0, 0) * pt.x + T.at<double>(0, 1) * pt.y + T.at<double>(0, 2);
res[1] = T.at<double>(1, 0) * pt.x + T.at<double>(1, 1) * pt.y + T.at<double>(1, 2);
res[2] = T.at<double>(2, 0) * pt.x + T.at<double>(2, 1) * pt.y + T.at<double>(2, 2);;

if (abs(res[2]) < 1e-12) {
throw std::runtime_error("Infinity point");
}
return {res[0]/res[2], res[1]/res[2]};
}

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