added some methods to the heightmap class to calculate steepness and statistics

This commit is contained in:
esther 2024-05-07 18:32:07 +02:00
parent ec0a31db5d
commit c47565eabb
3 changed files with 126 additions and 3 deletions

View File

@ -1,3 +1,4 @@
#include <cstddef>
#include <gdal/cpl_conv.h> #include <gdal/cpl_conv.h>
#include <iostream> #include <iostream>
#include <stdexcept> #include <stdexcept>
@ -5,7 +6,7 @@
#include <gdal/gdal.h> #include <gdal/gdal.h>
#include "gdal/gdal_priv.h" #include "gdal/gdal_priv.h"
#include <gdal/gdal_frmts.h> #include <gdal/gdal_frmts.h>
//#include "matplotlibcpp.h"
#include "contour_creator.hh" #include "contour_creator.hh"
@ -41,6 +42,7 @@ HeightMap::HeightMap(const char* filepath)
0, 0 ); 0, 0 );
if (error) { throw std::runtime_error("Could not read tif file!"); } if (error) { throw std::runtime_error("Could not read tif file!"); }
band->FlushCache(); band->FlushCache();
const int size = this->width*this->height;
} }
float HeightMap::get_pixel(int x, int y) float HeightMap::get_pixel(int x, int y)
{ {
@ -62,7 +64,7 @@ void HeightMap::blur(float standard_deviation)
*(kernel + i) = 1.0; *(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 parallel
{ {
#pragma omp for #pragma omp for
@ -91,3 +93,109 @@ void HeightMap::blur(float standard_deviation)
this->data = blurred; this->data = blurred;
blurred = nullptr; blurred = nullptr;
} }
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() {
}*/

View File

@ -20,6 +20,8 @@ class HeightMap
HeightMap(const char* filepath); HeightMap(const char* filepath);
float get_pixel(int x,int y); float get_pixel(int x,int y);
void blur(float standard_deviation); void blur(float standard_deviation);
void statistics();
void calculate_steepness();
}; };
/** /**

View File

@ -12,6 +12,8 @@
#include "gdal/gdal_priv.h" #include "gdal/gdal_priv.h"
#include <gdal/gdal_frmts.h> #include <gdal/gdal_frmts.h>
#include <omp.h> #include <omp.h>
//#include "matplotlibcpp.h"
//namespace plt = matplotlibcpp;
std::vector<Point> produce_cellmap(HeightMap* heightmap, float z) std::vector<Point> produce_cellmap(HeightMap* heightmap, float z)
{ {
@ -235,6 +237,17 @@ int main(int argc, const char* argv[])
map.blur(0.8); map.blur(0.8);
map.statistics();
map.calculate_steepness();
auto lines = create_lines(&map, 5); auto lines = create_lines(&map, 5);
write_output_file(lines, "out.geojson", &map); write_output_file(lines, "out.geojson", &map);
std::vector<int> y = {1, 2, 3};
//plt::plot(y, "bo-");
//plt::show();
return 0;
} }