From c47565eabbee25e5c1f26533554afbfa50baab9c Mon Sep 17 00:00:00 2001 From: esther Date: Tue, 7 May 2024 18:32:07 +0200 Subject: [PATCH 1/2] added some methods to the heightmap class to calculate steepness and statistics --- src/HeightMap.cpp | 114 +++++++++++++++++++++++++++++++++++++++-- src/contour_creator.hh | 2 + src/main.cpp | 13 +++++ 3 files changed, 126 insertions(+), 3 deletions(-) diff --git a/src/HeightMap.cpp b/src/HeightMap.cpp index 017a0a4..406e4cb 100644 --- a/src/HeightMap.cpp +++ b/src/HeightMap.cpp @@ -1,3 +1,4 @@ +#include #include #include #include @@ -5,7 +6,7 @@ #include #include "gdal/gdal_priv.h" #include - +//#include "matplotlibcpp.h" #include "contour_creator.hh" @@ -41,6 +42,7 @@ HeightMap::HeightMap(const char* filepath) 0, 0 ); if (error) { throw std::runtime_error("Could not read tif file!"); } band->FlushCache(); + const int size = this->width*this->height; } float HeightMap::get_pixel(int x, int y) { @@ -62,7 +64,7 @@ void HeightMap::blur(float standard_deviation) *(kernel + i) = 1.0; } - float* blurred = (float *) CPLMalloc(sizeof(float)*this->width*this->height); + float* blurred = new float[this->width*this->height]; #pragma omp parallel { #pragma omp for @@ -90,4 +92,110 @@ void HeightMap::blur(float standard_deviation) free(this->data); this->data = blurred; blurred = nullptr; -} \ No newline at end of file +} + +void HeightMap::statistics() +{ + + float avg = 0; + float var = 0; + float std = 0; + + for (int i = 0; i < this->width*this->height; i++) + {int x = i%this->width; + int y = i/this->width; + if (x<0 || x>=this->width) + { + x = 0; + } + + if (y<0 || y>=this->height) + { + y = 0; + } + avg += this->get_pixel(x, y);} + //std::cout << this->get_pixel(x, y);} + + avg = avg/(this->width*this->height); + std::cout << "Average value: " << avg <<"\n"; + + for (int i = 0; i < this->width*this->height; i++) + { + int x = i%this->width; + int y = i/this->width; + var += (avg - this->get_pixel(x, y))*(avg - this->get_pixel(x, y)); + } + std = sqrt(var)/(this->width*this->height-1); + var = var / (this->width*this->height-1); + std::cout << "std: " << std << "\n"; + std::cout << "var: " << var << "\n"; +}; + + +void HeightMap::calculate_steepness() +{ + int kernel_height = 5; + int kernel_width = 5; + int kernel_size = kernel_height*kernel_width; + + float* kernel = new float[sizeof(float)*kernel_size]; + float* steepness = new float[this->width*this->height]; + #pragma omp parallel + { + #pragma omp for + for (int i = 0; i < this->height * this->width; i++) + { + int x = i%this->width; + int y = i/this->width; + float max = 0; + float min = 0; + for (int j = 0; j < kernel_size; j++) + { + if (this->get_pixel(x, y) > max) + { + max = get_pixel(x, y); + } + + if (this->get_pixel(x, y) < min) + { + min = get_pixel(x, y); + } + + if (x<0 || x>=this->width) + { + x = 0; + } + + if (y<0 || y>=this->height) + { + y = 0; + } + + } + steepness[i] = max - min; + } + } + const char * filename = "steepness.tif"; + auto driver = GetGDALDriverManager()->GetDriverByName("GeoRaster"); + GDALDataset *file; + file = driver->Create("steepness.tif", this->width, this->height, 1, GDT_Float32, NULL); + GDALRegister_GTiff(); + + file->AddBand(GDT_Float32, NULL); + GDALRasterBand *band = file->GetRasterBand(0); + CPLErr error = band->RasterIO( GF_Write, 0, 0, this->width, this->height, + steepness, this->width, this->height, GDT_Float32, + 0, 0 ); + GDALClose(file); +} + + +//namespace plt = matplotlibcpp; + +//matplotlibcpp::plot(x, y); +//plt::loglog(x, y); +/* +int main() { + + +}*/ \ No newline at end of file diff --git a/src/contour_creator.hh b/src/contour_creator.hh index 84dc926..2fe808a 100644 --- a/src/contour_creator.hh +++ b/src/contour_creator.hh @@ -20,6 +20,8 @@ class HeightMap HeightMap(const char* filepath); float get_pixel(int x,int y); void blur(float standard_deviation); + void statistics(); + void calculate_steepness(); }; /** diff --git a/src/main.cpp b/src/main.cpp index 91e78e3..dfb7e09 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -12,6 +12,8 @@ #include "gdal/gdal_priv.h" #include #include +//#include "matplotlibcpp.h" +//namespace plt = matplotlibcpp; std::vector produce_cellmap(HeightMap* heightmap, float z) { @@ -234,7 +236,18 @@ int main(int argc, const char* argv[]) std::cout << "max: " << map.max << " min: " << map.min << "\n"; map.blur(0.8); + + map.statistics(); + + map.calculate_steepness(); auto lines = create_lines(&map, 5); write_output_file(lines, "out.geojson", &map); + + + std::vector y = {1, 2, 3}; + //plt::plot(y, "bo-"); + //plt::show(); + + return 0; } \ No newline at end of file From 97c52ecef3519fa22e1eeacb20596da8e5f5325e Mon Sep 17 00:00:00 2001 From: Trygve Date: Tue, 7 May 2024 19:38:34 +0200 Subject: [PATCH 2/2] Tried fixing gdal bullshit in steepnes method --- src/HeightMap.cpp | 20 ++++++-------------- src/contour_creator.hh | 2 +- 2 files changed, 7 insertions(+), 15 deletions(-) diff --git a/src/HeightMap.cpp b/src/HeightMap.cpp index 406e4cb..f3db3d3 100644 --- a/src/HeightMap.cpp +++ b/src/HeightMap.cpp @@ -33,6 +33,7 @@ HeightMap::HeightMap(const char* filepath) //https://gdal.org/api/gdaldataset_cpp.html#_CPPv4N11GDALDataset15GetGeoTransformEPd this->geotransform = (double *) CPLMalloc(sizeof(double)*6); this->reference_system = *(file->GetSpatialRef()); + this->filepath = filepath; file->GetGeoTransform(this->geotransform); @@ -175,11 +176,13 @@ void HeightMap::calculate_steepness() steepness[i] = max - min; } } + GDALAllRegister(); + GDALDataset *original_file = (GDALDataset *) GDALOpen( filepath, GA_ReadOnly ); + const char * filename = "steepness.tif"; auto driver = GetGDALDriverManager()->GetDriverByName("GeoRaster"); GDALDataset *file; - file = driver->Create("steepness.tif", this->width, this->height, 1, GDT_Float32, NULL); - GDALRegister_GTiff(); + file = driver->CreateCopy("steepness.tif", original_file, FALSE, NULL, NULL, NULL); file->AddBand(GDT_Float32, NULL); GDALRasterBand *band = file->GetRasterBand(0); @@ -187,15 +190,4 @@ void HeightMap::calculate_steepness() steepness, this->width, this->height, GDT_Float32, 0, 0 ); GDALClose(file); -} - - -//namespace plt = matplotlibcpp; - -//matplotlibcpp::plot(x, y); -//plt::loglog(x, y); -/* -int main() { - - -}*/ \ No newline at end of file +} \ No newline at end of file diff --git a/src/contour_creator.hh b/src/contour_creator.hh index 2fe808a..841ec0b 100644 --- a/src/contour_creator.hh +++ b/src/contour_creator.hh @@ -16,8 +16,8 @@ class HeightMap float min; //!< Minimum value in image float max; //!< Maximum value in image OGRSpatialReference reference_system; - HeightMap(const char* filepath); + const char* filepath; float get_pixel(int x,int y); void blur(float standard_deviation); void statistics();