Skip to content

Commit 56c5f24

Browse files
added gaussianWeightedAvg function
1 parent 28d9f75 commit 56c5f24

File tree

3 files changed

+68
-31
lines changed

3 files changed

+68
-31
lines changed

depthmap/depthmap.cpp

+62-22
Original file line numberDiff line numberDiff line change
@@ -1019,34 +1019,74 @@ void OrthoDepthmap::integratedCamera(const CameraDepthmap& camera, const char *o
10191019
}
10201020
//fit h = a+b*elev
10211021

1022-
void OrthoDepthmap::fitLinearRegressionFromPairs() {
1023-
float sumElevation = 0; // Somma delle elevazioni (x)
1024-
float sumH = 0; // Somma dei valori h (y)
1025-
float sumElevationH = 0; // Somma delle moltiplicazioni elevation * h
1026-
float sumElevationSquared = 0; // Somma delle elevazioni al quadrato
1027-
size_t n = point_cloud.size(); // Numero di punti nella nuvola di punti
1028-
1029-
// Ciclo sui punti per raccogliere i dati
1030-
for (const auto& point : point_cloud) {
1031-
float elevation = point[0]; // Elevation (x)
1032-
float h = point[2]; // Altezza (y)
1033-
1034-
sumElevation += elevation;
1035-
sumH += h;
1036-
sumElevationH += elevation * h;
1037-
sumElevationSquared += elevation * elevation;
1038-
}
1022+
void OrthoDepthmap::gaussianWeightedAvg(const char *textPath, int grid_x, int grid_y, float sigma) {
10391023

1040-
float b = (n * sumElevationH - sumElevation * sumH) / (n * sumElevationSquared - sumElevation * sumElevation);
1041-
float a = (sumH - b * sumElevation) / n;
1024+
loadPointCloud(textPath);
1025+
std::vector<float> x, y, z;
10421026

10431027
for (const auto& point : point_cloud) {
1044-
float elevation = point[0];
1045-
float hCalculated = a + b * elevation;
1046-
cout << "For elevation:" << elevation << "calculated h:" << hCalculated << endl;
1028+
x.push_back(point[0]);
1029+
y.push_back(point[1]);
1030+
z.push_back(point[2]);
1031+
}
1032+
1033+
float x_min = *std::min_element(x.begin(), x.end());
1034+
float x_max = *std::max_element(x.begin(), x.end());
1035+
float y_min = *std::min_element(y.begin(), y.end());
1036+
float y_max = *std::max_element(y.begin(), y.end());
1037+
1038+
float x_step = (x_max - x_min) / (grid_x - 1);
1039+
float y_step = (y_max - y_min) / (grid_y - 1);
1040+
1041+
std::vector<std::vector<float>> z_grid(grid_x, std::vector<float>(grid_y, std::numeric_limits<float>::quiet_NaN()));
1042+
float max_distance = 3 * sigma;
1043+
for (size_t p = 0; p < x.size(); p++) {
1044+
float px = x[p];
1045+
float py = y[p];
1046+
float pz = z[p];
1047+
int x_start = std::max(0, static_cast<int>((px - max_distance - x_min) / x_step));
1048+
int x_end = std::min(grid_x - 1, static_cast<int>((px + max_distance - x_min) / x_step));
1049+
int y_start = std::max(0, static_cast<int>((py - max_distance - y_min) / y_step));
1050+
int y_end = std::min(grid_y - 1, static_cast<int>((py + max_distance - y_min) / y_step));
1051+
1052+
for (int i = x_start; i <= x_end; i++) {
1053+
for (int j = y_start; j <= y_end; j++) {
1054+
float xg = x_min + i * x_step;
1055+
float yg = y_min + j * y_step;
1056+
1057+
float distance = std::sqrt((px - xg) * (px - xg) + (py - yg) * (py - yg));
1058+
if (distance <= max_distance) {
1059+
float weight = std::exp(-(distance * distance) / (2 * sigma * sigma));
1060+
1061+
if (qIsNaN(z_grid[i][j])) {
1062+
z_grid[i][j] = weight * pz;
1063+
} else {
1064+
z_grid[i][j] += weight * pz;
1065+
}
1066+
}
1067+
}
1068+
}
10471069
}
10481070

1071+
for (int i = 0; i < grid_x; i++) {
1072+
for (int j = 0; j < grid_y; j++) {
1073+
if (!std::isnan(z_grid[i][j])) {
1074+
z_grid[i][j] /= (3 * sigma);
1075+
}
1076+
}
1077+
}
1078+
for (int i = 0; i < grid_x; i++) {
1079+
for (int j = 0; j < grid_y; j++) {
1080+
if (!std::isnan(z_grid[i][j])) {
1081+
qDebug() << "z_grid[" << i << "][" << j << "] = " << z_grid[i][j];
1082+
} else {
1083+
qDebug() << "z_grid[" << i << "][" << j << "] = NaN";
1084+
}
1085+
}
1086+
}
10491087
}
1088+
1089+
10501090
//Real to Pixel Coordinates: (-0.981, 2.041, -11.132) -> Pixel: (110.688, 65.4375, -27.75)
10511091
//point outside the image limits Point 3D: (-0.981, 2.041, -11.132), Coordinate pixel: (-1, 2)
10521092

depthmap/depthmap.h

+1-6
Original file line numberDiff line numberDiff line change
@@ -90,12 +90,7 @@ class OrthoDepthmap:
9090
//legge nella depth l h corrispondente
9191
void verifyPointCloud();
9292
void integratedCamera(const CameraDepthmap& camera, const char *outputFile);
93-
void fitLinearRegressionFromPairs();
94-
95-
96-
97-
98-
93+
void gaussianWeightedAvg(const char *textPath, int grid_x, int grid_y, float sigma);
9994
/*1.
10095
*/
10196

depthmap/main.cpp

+5-3
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,7 @@ int main(int argc, char *argv[]) {
1414
return 1;
1515
}*/
1616
//input
17-
#define MACOS 1
17+
//#define MACOS 1
1818
#ifdef MACOS
1919
QString base = "/Users/erika/Desktop/";
2020
#else
@@ -25,7 +25,8 @@ int main(int argc, char *argv[]) {
2525
QString cameraDepthmap = base + "testcenterRel_copia/datasets/L04C12_depth_rti.tiff";
2626
QString orientationXmlPath = base + "testcenterRel_copia/photogrammetry/Ori-Relative/Orientation-L04C12.tif.xml";
2727
QString maskPath = base + "testcenterRel_copia/photogrammetry/Malt/Masq_STD-MALT_DeZoom4.tif";
28-
QString plyFile = base +"testcenterRel_copia/photogrammetry/AperiCloud_Relative__mini.ply";
28+
QString plyFile = base +"testcenterRel_copia/photogrammetry/AperiCloud_Relative_mini.ply";
29+
QString point_txt = base + "testcenterRel_copia/photogrammetry/points_h.txt";
2930
Depthmap depth;
3031

3132
//output
@@ -52,7 +53,6 @@ int main(int argc, char *argv[]) {
5253
exit(-1);
5354
}
5455
ortho.verifyPointCloud();
55-
ortho.fitLinearRegressionFromPairs();
5656

5757
CameraDepthmap depthCam;
5858
//xml per camera e tiff per la depth map
@@ -64,6 +64,8 @@ int main(int argc, char *argv[]) {
6464
}
6565
ortho.integratedCamera(depthCam, qPrintable(output_points));
6666

67+
ortho.gaussianWeightedAvg(qPrintable(point_txt), 40, 40, 0.025);
68+
6769

6870
//int pixelX = 165;
6971
//int pixelY = 144;

0 commit comments

Comments
 (0)