diff --git a/src/main.cpp b/src/main.cpp index 8da17d0..e3f9de3 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -36,55 +36,6 @@ std::vector produce_cellmap(HeightMap* heightmap, float z) } - - - -std::vector> vector_cellmap(HeightMap* heightmap, int interval) -{ - int num_contours = (heightmap->max - heightmap->min)/interval; - std::vector> vector_contours; - omp_set_num_threads(12); - - #pragma omp parallel - { - std::vector> vec_private; - #pragma omp for - 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++){ - Point current = points[j]; - //std::cout << points_allocated << " " << points.size() << "\n"; - for (int k = 0; k < points.size(); k++){ - Point candidate = points[k]; - if (candidate.next != nullptr) { - if (candidate.x +1 == current.x) { - current.next = &candidate; - } - else if (candidate.x -1 == current.x) { - current.next = &candidate; - } - else if (candidate.y +1 == current.x) { - current.next = &candidate; - } - else if (candidate.y -1 == current.x) { - current.next = &candidate; - } - - } - } - } - } - #pragma omp critical - - vector_contours.insert(vector_contours.end(), vec_private.begin(), vec_private.end()); - - } - return vector_contours; -} - std::tuple local_to_projected(double* geotransform, int x, int y) { return { @@ -94,6 +45,116 @@ std::tuple local_to_projected(double* geotransform, int x, int y geotransform[4] * 0.5 + geotransform[5] * 0.5 + geotransform[3] }; } +std::vector create_contours(OGRFeatureDefn* layer_def, std::vector points, HeightMap* heightmap) +{ + std::vector lines; + for (int i = 0; iSetField( "Name", i ); + Point* current = &points[0]; + int points_delegated = 0; + while (points_delegated < points.size()) + { + //int x_raw = i%width; + //int y_raw = i/height; + int x_raw = current->x; + int y_raw = current->y; + auto [x, y] = local_to_projected(heightmap->geotransform, x_raw, y_raw); + //std::cout << x << ", " << y << " "; + //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 ); + points_delegated++; + current = current->next; + if (current == nullptr) + { + for (int k = 0; k < points.size(); k++) + { + if (points[k].previous == nullptr) + { + current = &points[k]; + } + } + if (current == nullptr) + { + break; + } + + } + } + + if ( feature->SetGeometry(geometry) != OGRERR_NONE) + { + printf( "Failed to set geometry.\n" ); + exit( 1 ); + } + OGRGeometryFactory::destroyGeometry(geometry); + lines.push_back(feature); + } + +} + +std::vector vector_cellmap(OGRFeatureDefn* layer_def, HeightMap* heightmap, int interval) +{ + int num_contours = (heightmap->max - heightmap->min)/interval; + std::vector vector_contours; + omp_set_num_threads(12); + + #pragma omp parallel + { + std::vector lines; + #pragma omp for + 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++){ + Point* current = &points[j]; + //std::cout << points_allocated << " " << points.size() << "\n"; + for (int k = 0; k < points.size(); k++){ + Point* candidate = &points[k]; + if (candidate->previous == nullptr) { + + if (candidate->x +1 == current->x && candidate->y == current->y) { + current->next = candidate; + candidate->previous = current; + } + else if (candidate->x -1 == current->x && candidate->y == current->y) { + current->next = candidate; + candidate->previous = current; + } + else if (candidate->y +1 == current->y && candidate->x == current->x) { + current->next = candidate; + candidate->previous = current; + } + else if (candidate->y -1 == current->y && candidate->x == current->x) { + current->next = candidate; + candidate->previous = current; + } + + } + } + } + auto lines = create_lines(layer_def, points, heightmap); + + } + + #pragma omp critical + + vector_contours.insert(vector_contours.end(), lines.begin(), lines.end()); + + } + return vector_contours; +} + + constexpr std::tuple marching_squares_lookup(int mcase) { @@ -114,7 +175,9 @@ constexpr std::tuple marching_squares_lookup(int case 14: return {0, 0.5, 0.5, 0}; }; } -void write_output_file(std::vector> all_points, const char *filepath, HeightMap* heightmap) + + +void write_output_file(const char *filepath, HeightMap* heightmap) { const char *pszDriverName = "GeoJSON"; @@ -160,47 +223,19 @@ void write_output_file(std::vector> all_points, const char *f int height = heightmap->height -1; int size = width * height; - for (int j = 0; j < all_points.size(); j++) + // gets the lines from the heightmap. Runs the actual algorithm + auto lines = create_contours(poLayer->GetLayerDefn(), heightmap, 5); + + + for (int j = 0; j < lines.size(); 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++) - { - - //int x_raw = i%width; - //int y_raw = i/height; - int x_raw = points[i].x; - int y_raw = points[i].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 ) + if( poLayer->CreateFeature( lines[j] ) != OGRERR_NONE ) { printf( "Failed to create feature in shapefile.\n" ); exit( 1 ); } - OGRFeature::DestroyFeature( feature ); + OGRFeature::DestroyFeature( lines[j] ); } @@ -220,6 +255,6 @@ int main(int argc, const char* argv[]) auto cellmap = produce_cellmap(&map, 40); - auto cellmaps = vector_cellmap(&map, 5); - write_output_file(cellmaps, "out.geojson", &map); + + write_output_file("out.geojson", &map); } \ No newline at end of file