Skip to content

Commit

Permalink
start measuring accuracy
Browse files Browse the repository at this point in the history
  • Loading branch information
Lypho2012 committed May 17, 2024
1 parent 9113ce0 commit 8ea2d4d
Show file tree
Hide file tree
Showing 3 changed files with 44 additions and 85 deletions.
2 changes: 1 addition & 1 deletion opensfm/dense.py
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ def compute_depthmaps(
continue
mind, maxd = compute_depth_range(graph, reconstruction, shot, config)
arguments.append((data, neighbors[shot.id], mind, maxd, shot))
parallel_map(compute_depthmap_catched, arguments, 1)
parallel_map(compute_depthmap_catched, arguments, processes)
#for arg in arguments:
# compute_depthmap_catched(arg)

Expand Down
125 changes: 42 additions & 83 deletions opensfm/src/dense/src/depthmap.cc
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,6 @@
#include </home/czhang/OpenSfM/bsiCPP/bsi/BsiSigned.hpp>
#include </home/czhang/OpenSfM/bsiCPP/bsi/BsiUnsigned.hpp>

#include <cmath>

namespace dense {

static const double z_epsilon = 1e-8;
Expand Down Expand Up @@ -258,11 +256,10 @@ void DepthmapEstimator::ComputePatchMatchSample(

double PRECISION = 1000000;
for (int i = 0; i < patchmatch_iterations_; ++i) {
//std::cout << i << "\n";
//auto t31 = std::chrono::high_resolution_clock::now();
std::cout << i << "\n";
//std::cout << i << "\n";
PatchMatchForwardPass(result, true);
std::cout << i << " " << "forward pass: " << H_.size() << "\n";
//std::cout << i << " " << "forward pass: " << H_.size() << "\n";
//auto t32 = std::chrono::high_resolution_clock::now();
//duration = std::chrono::duration_cast<std::chrono::microseconds>(t32 - t31);
//std::cout << "PatchMatchForwardPass: " << duration.count() << "\n";
Expand All @@ -274,17 +271,17 @@ void DepthmapEstimator::ComputePatchMatchSample(
for (int k=0; k < H_.size(); k++) {
vec.emplace_back(static_cast<long>(H_[k](h, j)*PRECISION));

std::cout << H_[k](h, j) << " ";
//std::cout << H_[k](h, j) << " ";
}
std::cout <<"\n";
//std::cout <<"\n";
//std::cout << "vec size: " << vec.size() << "\n";
H_bsi.emplace_back(bsi.buildBsiAttributeFromVectorSigned(vec,0.5));
std::cout << "build bsi\n";
//std::cout << "build bsi\n";
vec.clear();
std::cout << "clear vec\n";
//std::cout << "clear vec\n";
}
}
std::cout << "i: ";
/*std::cout << "i: ";
for (int k=0; k < H_.size(); k++) {
std::cout << i_[k] << " ";
}
Expand All @@ -293,50 +290,20 @@ void DepthmapEstimator::ComputePatchMatchSample(
for (int k=0; k < H_.size(); k++) {
std::cout << j_[k] << " ";
}
std::cout << "\n";
std::cout << "\n";*/
std::cout << "H_ size: " << H_.size() << " " << i_.size() << " " << j_.size() << " " << H_bsi.size() << "\n";
BsiAttribute<uint64_t>* i_bsi = bsi.buildBsiAttributeFromVectorSigned(i_,0.5);
BsiAttribute<uint64_t>* j_bsi = bsi.buildBsiAttributeFromVectorSigned(j_,0.5);
/*BsiAttribute<uint64_t>* u_bsi = H_bsi[0]->multiplyBSI(j_bsi)->SUM(H_bsi[1]->multiplyBSI(i_bsi)->SUM(H_bsi[2]));
BsiAttribute<uint64_t>* v_bsi = H_bsi[3]->multiplyBSI(j_bsi)->SUM(H_bsi[4]->multiplyBSI(i_bsi)->SUM(H_bsi[5]));
BsiAttribute<uint64_t>* w_bsi = H_bsi[6]->multiplyBSI(j_bsi)->SUM(H_bsi[7]->multiplyBSI(i_bsi)->SUM(H_bsi[8]));*/
/*BsiAttribute<uint64_t>* u_bsi = H_bsi[0]->multiplyBSI(j_bsi);
std::cout << "u_bsi 1\n";
u_bsi = u_bsi->SUM(H_bsi[1]->multiplyBSI(i_bsi));
std::cout << "u_bsi 2\n";
u_bsi = u_bsi->SUM(H_bsi[2]);
std::cout << "u_bsi 3\n";
BsiAttribute<uint64_t>* v_bsi = H_bsi[3]->multiplyBSI(j_bsi);
std::cout << "v_bsi 1\n";
v_bsi = v_bsi->SUM(H_bsi[4]->multiplyBSI(i_bsi));
std::cout << "v_bsi 2\n";
v_bsi = v_bsi->SUM(H_bsi[5]);
std::cout << "v_bsi 3\n";
//std::cout << H_bsi[6]->rows << "\n";
BsiAttribute<uint64_t>* w_bsi = H_bsi[6]->multiplyBSI(j_bsi);
std::cout << "w_bsi 1\n";
w_bsi = w_bsi->SUM(H_bsi[7]->multiplyBSI(i_bsi));
std::cout << "w_bsi 2\n";
w_bsi = w_bsi->SUM(H_bsi[8]);
std::cout << "w_bsi 3\n";*/

auto t1 = std::chrono::high_resolution_clock::now();
/*BsiAttribute<uint64_t>* u_bsi = H_bsi[0]->multiplyWithBsiHorizontal(j_bsi)->sum_Horizontal(H_bsi[1]->multiplyWithBsiHorizontal(i_bsi)->sum_Horizontal(H_bsi[2]));
BsiAttribute<uint64_t>* v_bsi = H_bsi[3]->multiplyWithBsiHorizontal(j_bsi)->sum_Horizontal(H_bsi[4]->multiplyWithBsiHorizontal(i_bsi)->sum_Horizontal(H_bsi[5]));
BsiAttribute<uint64_t>* w_bsi = H_bsi[6]->multiplyWithBsiHorizontal(j_bsi)->sum_Horizontal(H_bsi[7]->multiplyWithBsiHorizontal(i_bsi)->sum_Horizontal(H_bsi[8]));*/

/*BsiAttribute<uint64_t>* u_bsi = H_bsi[0]->multiplication_Horizontal(j_bsi)->SUM(H_bsi[1]->multiplication_Horizontal(i_bsi)->SUM(H_bsi[2]));
BsiAttribute<uint64_t>* v_bsi = H_bsi[3]->multiplication_Horizontal(j_bsi)->SUM(H_bsi[4]->multiplication_Horizontal(i_bsi)->SUM(H_bsi[5]));
BsiAttribute<uint64_t>* w_bsi = H_bsi[6]->multiplication_Horizontal(j_bsi)->SUM(H_bsi[7]->multiplication_Horizontal(i_bsi)->SUM(H_bsi[8]));*/

BsiAttribute<uint64_t>* u_bsi = H_bsi[0]->multiplyWithBsiHorizontal(j_bsi)->SUM(H_bsi[1]->multiplyWithBsiHorizontal(i_bsi)->SUM(H_bsi[2]));
BsiAttribute<uint64_t>* v_bsi = H_bsi[3]->multiplyWithBsiHorizontal(j_bsi)->SUM(H_bsi[4]->multiplyWithBsiHorizontal(i_bsi)->SUM(H_bsi[5]));
BsiAttribute<uint64_t>* w_bsi = H_bsi[6]->multiplyWithBsiHorizontal(j_bsi)->SUM(H_bsi[7]->multiplyWithBsiHorizontal(i_bsi)->SUM(H_bsi[8]));

auto t2 = std::chrono::high_resolution_clock::now();
auto duration = std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1);
std::cout << "Calculate u,v,w: BSI: " << duration.count() << " Regular: " << uvw_time<< "\n";
std::cout << "U size: " << u_bsi->rows << "V size: " << v_bsi->rows << "W size: " << w_bsi->rows << "\n";
std::cout << "U size: " << u_bsi->rows << " V size: " << v_bsi->rows << " W size: " << w_bsi->rows << "\n";

double residuals_u = 0;
double mean_u = 0;
Expand All @@ -348,42 +315,34 @@ void DepthmapEstimator::ComputePatchMatchSample(
double mean_w = 0;
double squares_w = 0;
for (int k=0; k<H_.size(); k++) {
/*residuals_u += pow(u_bsi->getValue(k)/PRECISION - u_[k],2);
std::cout << k << "\n";
double residual_u = u_bsi->getValue(k)/PRECISION - u_[k];
residuals_u += residual_u*residual_u;
std::cout << "u: " << residuals_u << "\n";
mean_u += u_[k];
residuals_v += pow(v_bsi->getValue(k)/PRECISION - v_[k],2);
double residual_v = v_bsi->getValue(k)/PRECISION - v_[k];
residuals_v += residual_v*residual_v;
std::cout << "v: " << residuals_v << "\n";
mean_v += v_[k];
residuals_w += pow(w_bsi->getValue(k)/PRECISION - w_[k],2);
mean_w += w_[k];*/
std::cout << "u: " << u_bsi->getValue(k)/PRECISION << " " << u_[k] << " v: " << v_bsi->getValue(k)/PRECISION << " " << v_[k] << " w: " << w_bsi->getValue(k)/PRECISION << " " << w_[k] << "\n";
double residual_w = w_bsi->getValue(k)/PRECISION - w_[k];
residuals_w += residual_w*residual_w;
std::cout << "w: " << residuals_w << "\n";
mean_w += w_[k];
}
/*mean_u /= H_.size();
mean_u /= H_.size();
mean_v /= H_.size();
mean_w /= H_.size();
for (int k=0; k<H_.size(); k++) {
squares_u += pow(u_[k] - mean_u, 2);
squares_v += pow(v_[k] - mean_v, 2);
squares_w += pow(w_[k] - mean_w, 2);
double square_u = u_[k] - mean_u;
squares_u += square_u*square_u;
double square_v = v_[k] - mean_v;
squares_v += square_v*square_v;
double square_w = w_[k] - mean_w;
squares_w += square_w*square_w;
}
std::cout << "MSE u: " << residuals_u/H_.size() << " MSE v: " << residuals_v/H_.size() << " MSE w: " << residuals_w/H_.size() << "\n";
std::cout << "R2 u: " << 1 - residuals_u/H_.size()/squares_u << " R2 v: " << 1 - residuals_v/H_.size()/squares_v << " R2 w: " << 1 - residuals_w/H_.size()/squares_w << "\n";*/
std::cout << "R2 u: " << 1 - residuals_u/H_.size()/squares_u << " R2 v: " << 1 - residuals_v/H_.size()/squares_v << " R2 w: " << 1 - residuals_w/H_.size()/squares_w << "\n";

BsiAttribute<uint64_t>* u_bsi1 = H_bsi[0]->multiplyWithBsiHorizontal(j_bsi);
BsiAttribute<uint64_t>* u_bsi2 = H_bsi[1]->multiplyWithBsiHorizontal(i_bsi);
BsiAttribute<uint64_t>* u_bsi3 = u_bsi1->SUM(u_bsi2);

std::cout << "u1 size: " << u_bsi1->rows << " u2 size: "<< u_bsi2->rows << " u3 size: "<< u_bsi3->rows << "\n";
for (int k=0; k<H_.size(); k++) {
std::cout << "u1: " << u_bsi1->getValue(k)/PRECISION << " " << Hj_[k] << " u2: " << u_bsi2->getValue(k)/PRECISION << " " << Hi_[k] << " u3: " << u_bsi3->getValue(k)/PRECISION << " " << Hij_[k] << "\n";
}
/*for (int k=0; k<H_.size(); k++) {
std::cout << "u1: " << u_bsi1->getValue(k)/PRECISION << " " << Hj_[k] << "\n";
}
for (int k=0; k<H_.size(); k++) {
std::cout << " u2: " << u_bsi2->getValue(k)/PRECISION << " " << Hi_[k] << "\n";
}
for (int k=0; k<H_.size(); k++) {
std::cout << " u3: " << u_bsi3->getValue(k)/PRECISION << " " << Hij_[k] << "\n";
}*/
H_.clear();
u_.clear();
v_.clear();
Expand Down Expand Up @@ -483,7 +442,7 @@ void DepthmapEstimator::PatchMatchForwardPass(DepthmapEstimatorResult *result,
int hpz = (patch_size_ - 1) / 2;
for (int i = hpz; i < result->depth.rows - hpz; ++i) {
for (int j = hpz; j < result->depth.cols - hpz; ++j) {
if (H_.size() >= 1000000) return;
if (H_.size() >= 100) return;
PatchMatchUpdatePixel(result, i, j, adjacent, sample);
//std::cout << "Forward " << i << " " << j << " " << H_.size() << "\n";
}
Expand Down Expand Up @@ -643,25 +602,13 @@ float DepthmapEstimator::ComputePlaneImageScore(int i, int j,
const cv::Vec3f &plane,
int other) {
if (H_.size() >= 100) return 0;
auto t1 = std::chrono::high_resolution_clock::now();
cv::Matx33f H = PlaneInducedHomographyBaked(Kinvs_[0], Qs_[other], as_[other],
Ks_[other], plane);
int hpz = (patch_size_ - 1) / 2;
auto t1 = std::chrono::high_resolution_clock::now();
float u = H(0, 0) * j + H(0, 1) * i + H(0, 2);
float v = H(1, 0) * j + H(1, 1) * i + H(1, 2);
float w = H(2, 0) * j + H(2, 1) * i + H(2, 2);
auto t2 = std::chrono::high_resolution_clock::now();
auto duration = std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1);
uvw_time += duration.count();
H_.push_back(H);
u_.push_back(u);
v_.push_back(v);
w_.push_back(w);
i_.push_back(i);
j_.push_back(j);
Hj_.push_back(H(0, 0) * j);
Hi_.push_back(H(0, 1) * i);
Hij_.push_back(H(0, 0) * j + H(0, 1) * i);

if (w == 0.0) {
return -1.0f;
Expand All @@ -678,6 +625,18 @@ float DepthmapEstimator::ComputePlaneImageScore(int i, int j,
float im1_center = images_[0].at<unsigned char>(i, j);

NCCEstimator ncc;
auto t2 = std::chrono::high_resolution_clock::now();
auto duration = std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1);
uvw_time += duration.count();
H_.push_back(H);
u_.push_back(u);
v_.push_back(v);
w_.push_back(w);
i_.push_back(i);
j_.push_back(j);
Hj_.push_back(H(0, 0) * j);
Hi_.push_back(H(0, 1) * i);
Hij_.push_back(H(0, 0) * j + H(0, 1) * i);
for (int dy = -hpz; dy <= hpz; ++dy) {
for (int dx = -hpz; dx <= hpz; ++dx) {
float im1 = images_[0].at<unsigned char>(i + dy, j + dx);
Expand Down

0 comments on commit 8ea2d4d

Please sign in to comment.