mirror of
https://gitlab.com/Trygve/contour-creator.git
synced 2024-12-04 13:00:18 +00:00
Compare commits
13 Commits
95649f005b
...
9300fa729c
Author | SHA1 | Date | |
---|---|---|---|
9300fa729c | |||
62d62987d3 | |||
5845ae6693 | |||
47d0e29868 | |||
43f8b78830 | |||
08255875c1 | |||
ec0a31db5d | |||
|
32a0b18f8e | ||
|
ab96f9477e | ||
790cf44f8b | |||
df7cac1095 | |||
93fcbc6a39 | |||
1ec3e8c778 |
3
.gitmodules
vendored
Normal file
3
.gitmodules
vendored
Normal file
@ -0,0 +1,3 @@
|
||||
[submodule "extern/argh"]
|
||||
path = extern/argh
|
||||
url = https://github.com/adishavit/argh.git
|
@ -1,15 +1,22 @@
|
||||
cmake_minimum_required(VERSION 3.20)
|
||||
|
||||
project(
|
||||
contour-creator
|
||||
LANGUAGES CXX)
|
||||
|
||||
find_package(GDAL CONFIG REQUIRED)
|
||||
|
||||
add_executable(${PROJECT_NAME}
|
||||
src/HeightMap.cpp src/CellMap.cpp src/main.cpp
|
||||
)
|
||||
|
||||
# Argh is a simple argrument parser
|
||||
add_subdirectory(extern/argh)
|
||||
target_link_libraries(${PROJECT_NAME} PRIVATE argh)
|
||||
|
||||
# Gdal is used for geodata IO
|
||||
find_package(GDAL CONFIG REQUIRED)
|
||||
target_link_libraries(${PROJECT_NAME} PRIVATE GDAL::GDAL)
|
||||
|
||||
find_package(OpenMP)
|
||||
if(OpenMP_CXX_FOUND)
|
||||
target_link_libraries(${PROJECT_NAME} PUBLIC OpenMP::OpenMP_CXX)
|
||||
endif()
|
||||
|
||||
target_link_libraries(${PROJECT_NAME} GDAL::GDAL)
|
@ -3,4 +3,4 @@ rm -rf build
|
||||
mkdir build
|
||||
cmake -DCMAKE_BUILD_TYPE=Debug -B build
|
||||
cmake --build build --parallel
|
||||
build/contour-creator 'example_files/crop.tif'
|
||||
build/contour-creator --interval=1 --blur example_files/crop.tif
|
||||
|
BIN
documentation/kurver_ås.png
Normal file
BIN
documentation/kurver_ås.png
Normal file
Binary file not shown.
After Width: | Height: | Size: 1.5 MiB |
3092
documentation/ms_wikipedia.svg
Normal file
3092
documentation/ms_wikipedia.svg
Normal file
File diff suppressed because it is too large
Load Diff
After Width: | Height: | Size: 88 KiB |
@ -5,16 +5,31 @@ author:
|
||||
- Trygve og Esther
|
||||
|
||||
---
|
||||
# What are contours
|
||||
|
||||
# Output of our program
|
||||
pretty pictures
|
||||
# Output of our program:
|
||||
![Contour map of Ås](kurver_ås.png)
|
||||
|
||||
# Marching squares
|
||||
![By Nicoguaro; Adroitwhiz - File:Marching_squares_algorithm.svg, CC BY-SA 4.0, https://commons.wikimedia.org/w/index.php?curid=143418710](ms_wikipedia.svg)
|
||||
|
||||
# The logical flow of our program
|
||||
flowchart
|
||||
|
||||
# GDAL
|
||||
What is it?
|
||||
```cpp
|
||||
int main(int argc, const char* argv[])
|
||||
{
|
||||
const char* filepath = argv[1];
|
||||
HeightMap map(filepath);
|
||||
map.blur(0.8)
|
||||
auto lines = create_lines(&map, 5);
|
||||
write_output_file(cellmaps, "out.geojson", &map);
|
||||
}
|
||||
```
|
||||
|
||||
# Performance
|
||||
show the benchmarks
|
||||
| | Time |
|
||||
|--------------------|--------------------|
|
||||
| 12 threads | 41s |
|
||||
| 1 thread | 116s |
|
||||
|
||||
|
||||
# Problems we encountered
|
||||
|
||||
|
1
extern/argh
vendored
Submodule
1
extern/argh
vendored
Submodule
@ -0,0 +1 @@
|
||||
Subproject commit 431bf323acd3bb805483491ff277dac5d85772ad
|
@ -43,15 +43,10 @@ class Point
|
||||
int x;
|
||||
int y;
|
||||
int mscase;
|
||||
Point * next;
|
||||
Point * previous;
|
||||
bool endpoint = false;
|
||||
bool allocated = false;
|
||||
Point(int x, int y, int mscase){
|
||||
this -> x = x;
|
||||
this -> y = y;
|
||||
this -> mscase = mscase;
|
||||
this -> next = nullptr;
|
||||
this -> previous = nullptr;
|
||||
};
|
||||
};
|
||||
|
174
src/main.cpp
174
src/main.cpp
@ -1,20 +1,22 @@
|
||||
#include "contour_creator.hh"
|
||||
#include <algorithm>
|
||||
#include <array>
|
||||
#include <gdal/ogr_api.h>
|
||||
#include "ogrsf_frmts.h"
|
||||
#include <gdal/ogr_core.h>
|
||||
#include <gdal/ogr_feature.h>
|
||||
#include <gdal/ogr_geometry.h>
|
||||
#include <iostream>
|
||||
#include <iterator>
|
||||
#include <ostream>
|
||||
#include <string>
|
||||
#include <sys/types.h>
|
||||
#include <tuple>
|
||||
#include <vector>
|
||||
#include <cstdint>
|
||||
#include <gdal/gdal.h>
|
||||
#include "gdal/gdal_priv.h"
|
||||
#include <gdal/gdal_frmts.h>
|
||||
#include <omp.h>
|
||||
#include <fstream>
|
||||
#include "argh.h"
|
||||
|
||||
std::vector<Point> produce_cellmap(HeightMap* heightmap, float z)
|
||||
{
|
||||
@ -24,28 +26,31 @@ std::vector<Point> produce_cellmap(HeightMap* heightmap, float z)
|
||||
for (int i = 0; i<length; i++) {
|
||||
int y = i/(heightmap->width-1);
|
||||
int x = i%(heightmap->width-1);
|
||||
uint8_t result = (heightmap->get_pixel(x,y)>z) +
|
||||
(heightmap->get_pixel(x+1,y)>z)*2 +
|
||||
(heightmap->get_pixel(x+1,y+1)>z)*4 +
|
||||
(heightmap->get_pixel(x,y+1)>z)*8;
|
||||
uint8_t result = (heightmap->get_pixel(x,y)>z)*8 +
|
||||
(heightmap->get_pixel(x+1,y)>z)*4 +
|
||||
(heightmap->get_pixel(x+1,y+1)>z)*2 +
|
||||
(heightmap->get_pixel(x,y+1)>z);
|
||||
if (result != 0 and result != 15 ) {
|
||||
points.push_back(Point(x, y, result));
|
||||
|
||||
};
|
||||
}
|
||||
|
||||
return points;
|
||||
|
||||
|
||||
}
|
||||
|
||||
bool is_in(int value, std::array<int, 8> array)
|
||||
{
|
||||
return std::binary_search(array.begin(), array.end(), value);
|
||||
}
|
||||
|
||||
|
||||
std::vector<std::vector<Point>> vector_cellmap(HeightMap* heightmap, int interval)
|
||||
std::vector<std::vector<Point>> create_lines(HeightMap* heightmap, int interval)
|
||||
{
|
||||
int num_contours = (heightmap->max - heightmap->min)/interval;
|
||||
std::vector<std::vector<Point>> vector_contours;
|
||||
omp_set_num_threads(12);
|
||||
|
||||
std::array<int,8> north{4, 5,6,7,8,9,10,11};
|
||||
std::array<int,8> south{1,2,5,6,9,10,13,14};
|
||||
std::array<int,8> east{2,3,4,5,10,11,12,13};
|
||||
std::array<int,8> west{1,3,5,7,8,10,12,14};
|
||||
|
||||
#pragma omp parallel
|
||||
{
|
||||
@ -60,39 +65,44 @@ std::vector<std::vector<Point>> vector_cellmap(HeightMap* heightmap, int interva
|
||||
int points_allocated = 0;
|
||||
x = points[0].x;
|
||||
y = points[0].y;
|
||||
int current_case = points[0].mscase;
|
||||
points[0].allocated = true;
|
||||
line.push_back(points[0]);
|
||||
for (int j = 0; j < points.size(); j++){
|
||||
for (int k = 0; k < points.size(); k++){
|
||||
Point* candidate = &points[k];
|
||||
if (!candidate->allocated) {
|
||||
if (candidate->x == x +1 && candidate->y == y) {
|
||||
if (candidate->x == x +1 && candidate->y == y /*&& is_in(current_case, east)*/) {
|
||||
candidate->allocated = true;
|
||||
x = candidate->x;
|
||||
y = candidate->y;
|
||||
current_case = candidate->mscase;
|
||||
line.push_back(*candidate);
|
||||
points_allocated++;
|
||||
break;
|
||||
}
|
||||
else if (candidate->x == x -1 && candidate->y == y) {
|
||||
else if (candidate->x == x -1 && candidate->y == y /*&& is_in(current_case, west)*/) {
|
||||
x = candidate->x;
|
||||
y = candidate->y;
|
||||
current_case = candidate->mscase;
|
||||
candidate->allocated = true;
|
||||
line.push_back(*candidate);
|
||||
points_allocated++;
|
||||
break;
|
||||
}
|
||||
else if (candidate->x == x && candidate->y == y + 1) {
|
||||
else if (candidate->x == x && candidate->y == y + 1 /*&& is_in(current_case, north)*/) {
|
||||
x = candidate->x;
|
||||
y = candidate->y;
|
||||
current_case = candidate->mscase;
|
||||
candidate->allocated = true;
|
||||
line.push_back(*candidate);
|
||||
points_allocated++;
|
||||
break;
|
||||
}
|
||||
else if (candidate->x == x && candidate->y == y - 1) {
|
||||
else if (candidate->x == x && candidate->y == y - 1 /*&& is_in(current_case, south)*/) {
|
||||
x = candidate->x;
|
||||
y = candidate->y;
|
||||
current_case = candidate->mscase;
|
||||
candidate->allocated = true;
|
||||
line.push_back(*candidate);
|
||||
points_allocated++;
|
||||
@ -110,23 +120,20 @@ std::vector<std::vector<Point>> vector_cellmap(HeightMap* heightmap, int interva
|
||||
if (!candidate->allocated)
|
||||
{
|
||||
line.push_back(*candidate);
|
||||
candidate->allocated = true;
|
||||
x = candidate->x;
|
||||
y = candidate->y;
|
||||
current_case = candidate->mscase;
|
||||
points_allocated++;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
//vec_private.push_back(points);
|
||||
//vector_contours.push_back(points);
|
||||
//std::cout << vector_contours.size();
|
||||
}
|
||||
#pragma omp critical
|
||||
|
||||
vector_contours.insert(vector_contours.end(), vec_private.begin(), vec_private.end());
|
||||
|
||||
}
|
||||
return vector_contours;
|
||||
}
|
||||
@ -144,20 +151,35 @@ std::tuple<double, double> local_to_projected(double* geotransform, int x, int y
|
||||
constexpr std::tuple<double, double, double, double> marching_squares_lookup(int mcase)
|
||||
{
|
||||
switch (mcase) {
|
||||
case 1: return {0, 0.5, 0.5, 0};
|
||||
case 2: return {1, 0.5, 0.5, 0};
|
||||
case 3: return {0, 0.5, 1, 0};
|
||||
case 4: return {0.5, 1, 1, 0.5};
|
||||
case 5: return {0, 0, 0, 0}; //FIXME
|
||||
case 6: return {0.5, 1, 0.5, 0};
|
||||
case 7: return {0, 0.5, 0.5, 1};
|
||||
case 8: return {0, 0.5, 0.5, 1};
|
||||
case 9: return {0.5, 1, 0.5, 0};
|
||||
case 10: return {0, 0, 0, 0}; //FIXME
|
||||
case 11: return {0.5, 1, 1, 0.5};
|
||||
case 12: return {0, 0.5, 1, 0.5};
|
||||
case 13: return {0.5, 0, 1, 0.5};
|
||||
case 14: return {0, 0.5, 0.5, 0};
|
||||
case 1:
|
||||
return {0, 0.5, 0.5, 0};
|
||||
case 14:
|
||||
return {0.5, 0, 0, 0.5};
|
||||
case 2:
|
||||
return {0.5, 0, 1, 0.5};
|
||||
case 13:
|
||||
return {1, 0.5, 0.5, 0};
|
||||
case 3:
|
||||
return {0, 0.5, 1, 0.5};
|
||||
case 12:
|
||||
return {1, 0.5, 0, 0.5};
|
||||
case 4:
|
||||
return {0.5, 1, 1, 0.5};
|
||||
case 11:
|
||||
return {1, 0.5, 0.5, 1};
|
||||
case 5:
|
||||
case 10:
|
||||
return {0, 0, 0, 0}; //FIXME
|
||||
case 6:
|
||||
return {0.5, 1, 0.5, 0};
|
||||
case 9:
|
||||
return {0.5, 0, 0.5, 1};
|
||||
case 7:
|
||||
return {0, 0.5, 0.5, 1};
|
||||
case 8:
|
||||
return {0.5, 0.5, 0, 0.5};
|
||||
|
||||
|
||||
};
|
||||
}
|
||||
void write_output_file(std::vector<std::vector<Point>> all_points, const char *filepath, HeightMap* heightmap)
|
||||
@ -213,31 +235,42 @@ void write_output_file(std::vector<std::vector<Point>> all_points, const char *f
|
||||
feature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() );
|
||||
feature->SetField( "Name", j );
|
||||
|
||||
|
||||
std::vector<Point> points = all_points[j];
|
||||
|
||||
bool left_to_right = true;
|
||||
|
||||
if ((points[0].y - points[1].y) < 0)
|
||||
left_to_right = true;
|
||||
else
|
||||
left_to_right = false;
|
||||
|
||||
for (int k = 0; k < points.size(); k++)
|
||||
{
|
||||
//int x_raw = i%width;
|
||||
//int y_raw = i/height;
|
||||
bool left_to_right = true;
|
||||
|
||||
if ((points[k].x - points[k-1].x) < 0)
|
||||
left_to_right = false;
|
||||
bool top_to_bottom = true;
|
||||
if ((points[k].y - points[k-1].y) < 0)
|
||||
left_to_right = false;
|
||||
|
||||
|
||||
|
||||
int x_raw = points[k].x;
|
||||
int y_raw = points[k].y;
|
||||
/*
|
||||
if (x_raw > 500 || y_raw > 400)
|
||||
{
|
||||
std::cout << "🤯";std::cout.flush();
|
||||
continue;
|
||||
}
|
||||
*/
|
||||
auto [x, y] = local_to_projected(heightmap->geotransform, x_raw, y_raw);
|
||||
//auto [x_bl, y_bl, x_tr, y_tr] = marching_squares_lookup(*(cellmap->cells + i) );
|
||||
//std::cout << x_raw << ", " << y_raw << " ";std::cout.flush();
|
||||
//geometry->setPoint( geometry->getNumPoints(), x+x_tr*0.25, y+y_tr*0.25 );
|
||||
//geometry->setPoint( geometry->getNumPoints(), x+x_bl*0.25, y+y_bl*0.25 );
|
||||
geometry->setPoint(geometry->getNumPoints(),x ,y );
|
||||
|
||||
//current sometimes becomes random pointer?????
|
||||
auto [x1, y1, x2, y2] = marching_squares_lookup(points[k].mscase);
|
||||
|
||||
if (left_to_right or !top_to_bottom)
|
||||
{
|
||||
geometry->setPoint(geometry->getNumPoints(),x + (x1/6) ,y + (y1/6));
|
||||
geometry->setPoint(geometry->getNumPoints(),x + (x2/6) ,y + (y2/6));
|
||||
}
|
||||
else {
|
||||
geometry->setPoint(geometry->getNumPoints(),x + (x2/6) ,y + (y2/6));
|
||||
geometry->setPoint(geometry->getNumPoints(),x + (x1/6) ,y + (y1/6));
|
||||
}
|
||||
}
|
||||
|
||||
if ( feature->SetGeometry(geometry) != OGRERR_NONE)
|
||||
@ -253,22 +286,39 @@ void write_output_file(std::vector<std::vector<Point>> all_points, const char *f
|
||||
}
|
||||
OGRFeature::DestroyFeature( feature );
|
||||
}
|
||||
|
||||
GDALClose( poDS );
|
||||
}
|
||||
|
||||
int main(int argc, const char* argv[])
|
||||
{
|
||||
const char* filepath = argv[1];
|
||||
HeightMap map(filepath);
|
||||
argh::parser cmdl(argv);
|
||||
|
||||
std::cout << "x: " << map.width << " y: " << map.height << "\n";
|
||||
std::cout << "max: " << map.max << " min: " << map.min << "\n";
|
||||
if (cmdl[{ "-h", "--help" }])
|
||||
{
|
||||
std::cout << "Usage:\n"
|
||||
<< "contour_creator [OPTIONS] <input_file>\n"
|
||||
<< "-o; --output <FILENAME.geojson> - File to write output to (Default: contours.geojson)\n"
|
||||
<< "-i; --interval <int> - Set the interval between contours (Default: 5)\n"
|
||||
<< "-b; --blur - Blur the image\n"
|
||||
<< "--stats - Print statistical information about the heightmap\n";
|
||||
exit(0);
|
||||
}
|
||||
int interval;
|
||||
cmdl({"-i", "--interval"}, 5) >> interval;
|
||||
if (interval <= 0)
|
||||
std::cerr << "Interval must be valid positive integer!" << "\n";
|
||||
|
||||
map.blur(0.8);
|
||||
std::string output_file;
|
||||
cmdl({"-o", "--output"}, "contours.geojson") >> output_file;
|
||||
|
||||
auto cellmap = produce_cellmap(&map, 40);
|
||||
std::string input_file;
|
||||
cmdl(1) >> input_file;
|
||||
|
||||
auto cellmaps = vector_cellmap(&map, 5);
|
||||
write_output_file(cellmaps, "out.geojson", &map);
|
||||
|
||||
HeightMap map(input_file.c_str());
|
||||
if (cmdl[{"-b", "--blur"}])
|
||||
map.blur(0.8);
|
||||
auto lines = create_lines(&map, interval);
|
||||
write_output_file(lines, output_file.c_str(), &map);
|
||||
std::cout << "Contours written to " << output_file << " 🗺️\n";
|
||||
}
|
Loading…
Reference in New Issue
Block a user