Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
90 changes: 58 additions & 32 deletions src/phg/mvs/depth_maps/pm_depth_maps.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -46,18 +46,20 @@ vector3d project(const vector3d& global_point, const phg::Calibration& calibrati
return pixel_with_depth;
}

// TODO 101 реализуйте unproject (вам поможет тест на идемпотентность project -> unproject в test_depth_maps_pm)
vector3d unproject(const vector3d& pixel, const phg::Calibration& calibration, const matrix34d& PtoWorld)
vector3d unproject(const vector3d& pix, const phg::Calibration& calib, const matrix34d& worldFromPose)
{
double depth = pixel[2]; // на самом деле это не глубина, это координата по оси +Z (вдоль которой смотрит камера в ее локальной системе координат)
double z_cam = pix[2]; // локальная Z-координата (вдоль оптической оси камеры)

vector3d local_point; // TODO 102 пустите луч pixel из calibration а затем возьмите ан нем точку у которой по оси +Z координата=depth
vector3d tmp = calib.unproject({ pix[0], pix[1] });

vector3d global_point; // TODO 103 переведите точку из локальной системы в глобальную
const double factor = z_cam / tmp[2];
vector3d local = tmp * factor;

return global_point;
vector3d global = worldFromPose * homogenize(local);
return global;
}


void PMDepthMapsBuilder::buildDepthMap(unsigned int camera_key, cv::Mat& depth_map_res, cv::Mat& normal_map_res, cv::Mat& cost_map_res, float depth_min, float depth_max)
{
rassert(camera_key < ncameras, 238192841294108);
Expand Down Expand Up @@ -112,16 +114,18 @@ void PMDepthMapsBuilder::refinement()
n0 = normal_map.at<vector3f>(j, i);

// 2) случайной пертурбации текущей гипотезы (мутация и уточнение того что уже смогли найти)
dp = r.nextf(d0 * 0.5f, d0 * 1.5); // TODO 104: сделайте так чтобы отклонение было тем меньше, чем номер итерации ближе к NITERATIONS, улучшило ли это результат?
np = cv::normalize(n0 + randomNormalObservedFromCamera(cameras_RtoWorld[ref_cam], r) * 0.5); // TODO 105: сделайте так чтобы отклонение было тем меньше, чем номер итерации ближе к NITERATIONS, улучшило ли это результат?
float t = float(iter) / NITERATIONS; // t = 0..1, где 1 — финальная итерация
float min_factor = std::max(0.1f, 0.5f * (1.0f - t));
float max_factor = std::max(0.2f, 1.5f * (1.0f - t));
dp = r.nextf(d0 * min_factor, d0 * max_factor);

np = cv::normalize(n0 + randomNormalObservedFromCamera(cameras_RtoWorld[ref_cam], r) * 0.5);

dp = std::max(ref_depth_min, std::min(ref_depth_max, dp));

// 3) новой случайной гипотезы из фрустума поиска (новые идеи, вечный поиск во всем пространстве)
// TODO 106: создайте случайную гипотезу dr+nr, вам поможет:
// - r.nextf(...)
// - ref_depth_min, ref_depth_max
// - randomNormalObservedFromCamera - поможет создать нормаль которая гарантированно смотрит на нас
dr = r.nextf(ref_depth_min, ref_depth_max);
nr = randomNormalObservedFromCamera(cameras_RtoWorld[ref_cam], r);
}

float best_depth = d0;
Expand Down Expand Up @@ -315,10 +319,13 @@ float PMDepthMapsBuilder::estimateCost(ptrdiff_t i, ptrdiff_t j, double d, const
vector3d pixel(i + 0.5, j + 0.5, d);
vector3d global_point = unproject(pixel, calibration, cameras_PtoWorld[ref_cam]);

if (!(i - COST_PATCH_RADIUS >= 0 && i + COST_PATCH_RADIUS < width))
if (!(i - COST_PATCH_RADIUS >= 0 && i + COST_PATCH_RADIUS < width)) {
return NO_COST;
if (!(j - COST_PATCH_RADIUS >= 0 && j + COST_PATCH_RADIUS < height))
}

if (!(j - COST_PATCH_RADIUS >= 0 && j + COST_PATCH_RADIUS < height)) {
return NO_COST;
}

std::vector<float> patch0, patch1;

Expand Down Expand Up @@ -348,35 +355,48 @@ float PMDepthMapsBuilder::estimateCost(ptrdiff_t i, ptrdiff_t j, double d, const
double x = neighb_proj[0];
double y = neighb_proj[1];

// TODO 205: замените этот наивный вариант nearest neighbor сэмплирования текстуры на билинейную интерполяцию (учтите что центр пикселя - .5 после запятой)
ptrdiff_t u = x;
ptrdiff_t v = y;
ptrdiff_t u = static_cast<ptrdiff_t>(x);
ptrdiff_t v = static_cast<ptrdiff_t>(y);

// TODO 108: добавьте проверку "попали ли мы в камеру номер neighb_cam?" если не попали - возвращаем NO_COST
if (u < 0 || u >= width || v < 0 || v >= height) {
return NO_COST;
}

float intensity = cameras_imgs_grey[neighb_cam].at<unsigned char>(v, u) / 255.0f;
patch1.push_back(intensity);
}
}

// TODO 109: реализуйте ZNCC https://en.wikipedia.org/wiki/Cross-correlation#Zero-normalized_cross-correlation_(ZNCC)
// или слайд #25 в лекции 5 про SGM и Cost-функции - https://my.compscicenter.ru/attachments/classes/slides_w2n8WNLY/photogrammetry_lecture_090321.pdf
// ZNCC
rassert(patch0.size() == patch1.size(), 12489185129326);
size_t n = patch0.size();
float mean0 = 0.0f;
float mean1 = 0.0f;
// ...

for (size_t k = 0; k < n; ++k) {
float a = patch0[k];
float b = patch1[k];
mean0 += a;
mean1 += b;
// ...
mean0 += patch0[k];
mean1 += patch1[k];
}
mean0 /= n;
mean1 /= n;
// ...

float sum_prod = 0.0f;
float sum_sq0 = 0.0f;
float sum_sq1 = 0.0f;

for (size_t k = 0; k < n; ++k) {
float a = patch0[k] - mean0;
float b = patch1[k] - mean1;
sum_prod += a * b;
sum_sq0 += a * a;
sum_sq1 += b * b;
}

float zncc = 0.0f;
float denom = std::sqrt(sum_sq0 * sum_sq1);
if (denom > 1e-6f) {
zncc = sum_prod / denom;
}

// ZNCC в диапазоне [-1; 1], 1: идеальное совпадение, -1: ничего общего
rassert(zncc == zncc, 23141241210380); // проверяем что не nan
Expand All @@ -393,8 +413,9 @@ float PMDepthMapsBuilder::estimateCost(ptrdiff_t i, ptrdiff_t j, double d, const

float PMDepthMapsBuilder::avgCost(std::vector<float>& costs)
{
if (costs.size() == 0)
if (costs.size() == 0) {
return NO_COST;
}

std::sort(costs.begin(), costs.end());

Expand All @@ -403,10 +424,15 @@ float PMDepthMapsBuilder::avgCost(std::vector<float>& costs)
float cost_sum = best_cost;
float cost_w = 1.0f;

// TODO 110 реализуйте какое-то "усреднение cost-ов по всем соседям", с ограничением что участвуют только COSTS_BEST_K_LIMIT лучших
// TODO 111 добавьте к этому усреднению еще одно ограничение: если cost больше чем best_cost*COSTS_K_RATIO - то такой cost подозрительно плохой и мы его не хотим учитывать (вероятно occlusion)
// TODO 112 а что если в пикселе occlusion, но best_cost - большой и поэтому отсечение по best_cost*COSTS_K_RATIO не срабатывает? можно ли это отсечение как-то выправить для такого случая?
// TODO 207 а что если добавить какой-нибудь бонус в случае если больше чем Х камер засчиталось? улучшается/ухудшается ли от этого что-то на herzjezu25? а при большем числе фотографий
const int k_limit = std::min((int)costs.size(), COSTS_BEST_K_LIMIT);
for (int i = 1; i < k_limit; ++i) {
float c = costs[i];
if (c > best_cost * COSTS_K_RATIO) {
break; // последующие cost'ы только больше (так как массив отсортирован)
}
cost_sum += c;
cost_w += 1.0f;
}

float avg_cost = cost_sum / cost_w;
return avg_cost;
Expand Down
68 changes: 34 additions & 34 deletions tests/test_depth_maps_pm.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -34,43 +34,43 @@
TEST(test_depth_maps_pm, SingleDepthMap)
{
// TODO этот код надо раскомментировать чтобы запустить тестирование:
// Dataset dataset = loadDataset(DATASET_DIR, DATASET_DOWNSCALE);
// phg::PMDepthMapsBuilder builder(dataset.ncameras, dataset.cameras_imgs, dataset.cameras_imgs_grey, dataset.cameras_labels, dataset.cameras_P, dataset.calibration);
//
// size_t ci = 2; // строим карту глубины для третьей камеры (индексация с нуля)
// size_t cameras_limit = 5; // учитывая первые пять фотографий датасета, т.е. две камеры слева и две камеры справа
//
// dataset.ncameras = cameras_limit;
// cv::Mat depth_map, normal_map, cost_map;
// builder.buildDepthMap(ci, depth_map, cost_map, normal_map, dataset.cameras_depth_min[ci], dataset.cameras_depth_max[ci]);
Dataset dataset = loadDataset(DATASET_DIR, DATASET_DOWNSCALE);
phg::PMDepthMapsBuilder builder(dataset.ncameras, dataset.cameras_imgs, dataset.cameras_imgs_grey, dataset.cameras_labels, dataset.cameras_P, dataset.calibration);

size_t ci = 2; // строим карту глубины для третьей камеры (индексация с нуля)
size_t cameras_limit = 5; // учитывая первые пять фотографий датасета, т.е. две камеры слева и две камеры справа

dataset.ncameras = cameras_limit;
cv::Mat depth_map, normal_map, cost_map;
builder.buildDepthMap(ci, depth_map, cost_map, normal_map, dataset.cameras_depth_min[ci], dataset.cameras_depth_max[ci]);
}

TEST(test_depth_maps_pm, AllDepthMaps)
{
// TODO этот код можно раскомментировать чтобы построить много карт глубины и сохранить их облака точек:
// Dataset full_dataset = loadDataset(DATASET_DIR, DATASET_DOWNSCALE);
//
// const size_t ref_camera_shift = 2;
// const size_t to_shift = 5;
//
// std::vector<cv::Vec3d> all_points;
// std::vector<cv::Vec3b> all_colors;
// std::vector<cv::Vec3d> all_normals;
//
// size_t ndepth_maps = 0;
//
// for (size_t from = 0; from + to_shift <= full_dataset.ncameras; ++from) {
// size_t to = from + to_shift;
//
// Dataset dataset = full_dataset.subset(from, to);
//
// phg::PMDepthMapsBuilder builder(dataset.ncameras, dataset.cameras_imgs, dataset.cameras_imgs_grey, dataset.cameras_labels, dataset.cameras_P, dataset.calibration);
// cv::Mat depth_map, normal_map, cost_map;
// builder.buildDepthMap(ref_camera_shift, depth_map, cost_map, normal_map, dataset.cameras_depth_min[ref_camera_shift], dataset.cameras_depth_max[ref_camera_shift]);
// phg::PMDepthMapsBuilder::buildGoodPoints(depth_map, normal_map, cost_map, dataset.cameras_imgs[ref_camera_shift], dataset.calibration, builder.getCameraPtoWorld(ref_camera_shift), all_points, all_colors, all_normals);
// ++ndepth_maps;
//
// std::string tie_points_filename = std::string("data/debug/test_depth_maps_pm/") + getTestName() + "/all_points_" + to_string(ndepth_maps) + ".ply";
// phg::exportPointCloud(all_points, tie_points_filename, all_colors, all_normals);
// }
Dataset full_dataset = loadDataset(DATASET_DIR, DATASET_DOWNSCALE);

const size_t ref_camera_shift = 2;
const size_t to_shift = 5;

std::vector<cv::Vec3d> all_points;
std::vector<cv::Vec3b> all_colors;
std::vector<cv::Vec3d> all_normals;

size_t ndepth_maps = 0;

for (size_t from = 0; from + to_shift <= full_dataset.ncameras; ++from) {
size_t to = from + to_shift;

Dataset dataset = full_dataset.subset(from, to);

phg::PMDepthMapsBuilder builder(dataset.ncameras, dataset.cameras_imgs, dataset.cameras_imgs_grey, dataset.cameras_labels, dataset.cameras_P, dataset.calibration);
cv::Mat depth_map, normal_map, cost_map;
builder.buildDepthMap(ref_camera_shift, depth_map, cost_map, normal_map, dataset.cameras_depth_min[ref_camera_shift], dataset.cameras_depth_max[ref_camera_shift]);
phg::PMDepthMapsBuilder::buildGoodPoints(depth_map, normal_map, cost_map, dataset.cameras_imgs[ref_camera_shift], dataset.calibration, builder.getCameraPtoWorld(ref_camera_shift), all_points, all_colors, all_normals);
++ndepth_maps;

std::string tie_points_filename = std::string("data/debug/test_depth_maps_pm/") + getTestName() + "/all_points_" + to_string(ndepth_maps) + ".ply";
phg::exportPointCloud(all_points, tie_points_filename, all_colors, all_normals);
}
}
Loading