diff --git a/documentation/kurver_ås.png b/documentation/kurver_ås.png new file mode 100644 index 0000000..337cbfd Binary files /dev/null and b/documentation/kurver_ås.png differ diff --git a/documentation/ms_wikipedia.svg b/documentation/ms_wikipedia.svg new file mode 100644 index 0000000..f0a0f15 --- /dev/null +++ b/documentation/ms_wikipedia.svg @@ -0,0 +1,3092 @@ + + + + + + + 1 + 2 + 3 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 3 + 3 + 3 + 3 + 2 + 2 + 2 + + + + + + + + + + + + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 0 + 0 + 0 + 0 + 0 + 0 + 0 + 0 + 0 + 0 + 0 + 0 + 0 + 0 + 0 + 0 + + + + + + + + + + + + + + 1 + 2 + 3 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 3 + 3 + 3 + 3 + 2 + 2 + 2 + + + + + + + + + + + + Case 0 + Case 1 + Case 2 + Case 3 + Case 4 + Case 5 + Case 6 + Case 7 + Case 8 + Case 9 + Case 10 + Case 11 + Case 12 + Case 13 + Case 14 + Case 15 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1 + 2 + 4 + 8 + + + Give every cell a bitnumber based onwhich corners aretrue/false + Look up the contourlines in the tableand put them inthe cells + Look at the originalvalues and use linearinterpolation todetermine amore accurate positionof all the line end-points + Thresholdwith isovalue + here:isovalue = 1 + 0 = below, 1 = above + = below + =above + Binary imageto cells + Look-up table contour lines + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + 3 + 3 + 1 + 6 + 15 + 15 + 9 + 6 + 4 + 15 + 12 + 15 + 12 + 8 + 9 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/documentation/slides.md b/documentation/slides.md index a91d102..eacce3f 100644 --- a/documentation/slides.md +++ b/documentation/slides.md @@ -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 \ No newline at end of file +| | Time | +|--------------------|--------------------| +| 12 threads | 41s | +| 1 thread | 116s | + + +# Problems we encountered + diff --git a/src/contour_creator.hh b/src/contour_creator.hh index 4bc45f5..84dc926 100644 --- a/src/contour_creator.hh +++ b/src/contour_creator.hh @@ -43,13 +43,10 @@ class Point int x; int y; int mscase; - Point * next; - Point * previous; + 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; }; }; diff --git a/src/main.cpp b/src/main.cpp index 511f163..91e78e3 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -5,7 +5,6 @@ #include #include #include -#include #include #include #include @@ -22,28 +21,21 @@ std::vector produce_cellmap(HeightMap* heightmap, float z) for (int i = 0; iwidth-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; - - } - - -std::vector> vector_cellmap(HeightMap* heightmap, int interval) +std::vector> create_lines(HeightMap* heightmap, int interval) { int num_contours = (heightmap->max - heightmap->min)/interval; - std::vector> vector_contours; // initiating a vector of contours - omp_set_num_threads(12); + std::vector> vector_contours; #pragma omp parallel { @@ -52,26 +44,67 @@ std::vector> vector_cellmap(HeightMap* heightmap, int interva for (int i = 1; i <= num_contours; i++) { auto points = produce_cellmap(heightmap, heightmap->min + interval*i); - - for (int j = 0; j < points.size(); j++){ //looping through all the points. - Point current = points[j]; // assigning an instance of Point class to a point. - //std::cout << points_allocated << " " << points.size() << "\n"; - for (int k = 0; k < points.size(); k++){ // looping through all other points. - Point candidate = points[k]; - if (candidate.next != nullptr) { - if (candidate.x +1 == current.x) { - current.next = &candidate; + std::vector line; + int x; + int y; + int points_allocated = 0; + x = points[0].x; + y = points[0].y; + 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) { + candidate->allocated = true; + x = candidate->x; + y = candidate->y; + line.push_back(*candidate); + points_allocated++; + break; } - else if (candidate.x -1 == current.x) { - current.next = &candidate; + else if (candidate->x == x -1 && candidate->y == y) { + x = candidate->x; + y = candidate->y; + candidate->allocated = true; + line.push_back(*candidate); + points_allocated++; + break; } - else if (candidate.y +1 == current.x) { - current.next = &candidate; + else if (candidate->x == x && candidate->y == y + 1) { + x = candidate->x; + y = candidate->y; + candidate->allocated = true; + line.push_back(*candidate); + points_allocated++; + break; } - else if (candidate.y -1 == current.x) { - current.next = &candidate; + else if (candidate->x == x && candidate->y == y - 1) { + x = candidate->x; + y = candidate->y; + candidate->allocated = true; + line.push_back(*candidate); + points_allocated++; + break; + } + } + } + if (j > points_allocated) + { + vec_private.push_back(line); + line.clear(); + //std::cout << points.size() << " "; + for (int k = 0; k < points.size(); k++){ + Point* candidate = &points[k]; + if (!candidate->allocated) + { + line.push_back(*candidate); + x = candidate->x; + y = candidate->y; + points_allocated++; + break; } - } } } @@ -98,7 +131,7 @@ constexpr std::tuple marching_squares_lookup(int 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 3: return {0, 0.5, 1, 0.5}; 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}; @@ -138,7 +171,7 @@ void write_output_file(std::vector> all_points, const char *f OGRLayer *poLayer; - poLayer = poDS->CreateLayer( "contours", /*&(cellmaps[0]).reference_system*/ NULL, wkbLineString, NULL ); + poLayer = poDS->CreateLayer( "contours", &heightmap->reference_system, wkbLineString, NULL ); if( poLayer == NULL ) { printf( "Layer creation failed.\n" ); @@ -160,49 +193,35 @@ void write_output_file(std::vector> all_points, const char *f for (int j = 0; j < all_points.size(); j++) { - OGRFeature *feature; - OGRLineString *geometry = new OGRLineString(); - feature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() ); - feature->SetField( "Name", j ); - + OGRFeature *feature; + OGRLineString *geometry = new OGRLineString(); + feature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() ); + feature->SetField( "Name", j ); std::vector points = all_points[j]; - for (int i = 0; i < points.size(); i++) + for (int k = 0; k < points.size(); k++) { - - //int x_raw = i%width; - //int y_raw = i/height; - int x_raw = points[i].x; - int y_raw = points[i].y; + int x_raw = points[k].x; + int y_raw = points[k].y; 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 << ", " << y << " "; - //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 ); - - - - } - if ( feature->SetGeometry(geometry) != OGRERR_NONE) - { - printf( "Failed to set geometry.\n" ); - exit( 1 ); - } - OGRGeometryFactory::destroyGeometry(geometry); - if( poLayer->CreateFeature( feature ) != OGRERR_NONE ) - { - printf( "Failed to create feature in shapefile.\n" ); - exit( 1 ); - } - OGRFeature::DestroyFeature( feature ); - - + if ( feature->SetGeometry(geometry) != OGRERR_NONE) + { + printf( "Failed to set geometry.\n" ); + exit( 1 ); + } + OGRGeometryFactory::destroyGeometry(geometry); + if( poLayer->CreateFeature( feature ) != OGRERR_NONE ) + { + printf( "Failed to create feature in shapefile.\n" ); + exit( 1 ); + } + OGRFeature::DestroyFeature( feature ); } - GDALClose( poDS ); } @@ -216,8 +235,6 @@ int main(int argc, const char* argv[]) map.blur(0.8); - auto cellmap = produce_cellmap(&map, 40); - - auto cellmaps = vector_cellmap(&map, 5); - write_output_file(cellmaps, "out.geojson", &map); + auto lines = create_lines(&map, 5); + write_output_file(lines, "out.geojson", &map); } \ No newline at end of file