mirror of
				https://gitlab.com/Trygve/contour-creator.git
				synced 2025-11-04 02:40:45 +00:00 
			
		
		
		
	Seperate line creation so it can be paralized
This commit is contained in:
		
							parent
							
								
									77467594bd
								
							
						
					
					
						commit
						e07ab1445a
					
				
							
								
								
									
										205
									
								
								src/main.cpp
									
									
									
									
									
								
							
							
						
						
									
										205
									
								
								src/main.cpp
									
									
									
									
									
								
							@ -36,55 +36,6 @@ std::vector<Point> produce_cellmap(HeightMap* heightmap, float z)
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
std::vector<std::vector<Point>> vector_cellmap(HeightMap* heightmap, int interval)
 | 
			
		||||
{
 | 
			
		||||
    int num_contours = (heightmap->max - heightmap->min)/interval;
 | 
			
		||||
    std::vector<std::vector<Point>> vector_contours;
 | 
			
		||||
    omp_set_num_threads(12);
 | 
			
		||||
 | 
			
		||||
    #pragma omp parallel
 | 
			
		||||
    {
 | 
			
		||||
        std::vector<std::vector<Point>> 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<double, double> local_to_projected(double* geotransform, int x, int y)
 | 
			
		||||
{
 | 
			
		||||
    return {
 | 
			
		||||
@ -94,6 +45,116 @@ std::tuple<double, double> local_to_projected(double* geotransform, int x, int y
 | 
			
		||||
        geotransform[4] * 0.5 + geotransform[5] * 0.5 + geotransform[3]
 | 
			
		||||
        };
 | 
			
		||||
}
 | 
			
		||||
std::vector<OGRFeature*> create_contours(OGRFeatureDefn* layer_def, std::vector<Point> points, HeightMap* heightmap)
 | 
			
		||||
{
 | 
			
		||||
    std::vector<OGRFeature*> lines;
 | 
			
		||||
    for (int i = 0; i<points.size();i++)
 | 
			
		||||
    {
 | 
			
		||||
        OGRFeature *feature;
 | 
			
		||||
            OGRLineString *geometry = new OGRLineString();
 | 
			
		||||
            
 | 
			
		||||
            feature = OGRFeature::CreateFeature(layer_def);
 | 
			
		||||
            feature->SetField( "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<OGRFeature*> vector_cellmap(OGRFeatureDefn* layer_def, HeightMap* heightmap, int interval)
 | 
			
		||||
{
 | 
			
		||||
    int num_contours = (heightmap->max - heightmap->min)/interval;
 | 
			
		||||
    std::vector<OGRFeature*> vector_contours;
 | 
			
		||||
    omp_set_num_threads(12);
 | 
			
		||||
 | 
			
		||||
    #pragma omp parallel
 | 
			
		||||
    {
 | 
			
		||||
        std::vector<OGRFeature*> 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<double, double, double, double> marching_squares_lookup(int mcase)
 | 
			
		||||
{
 | 
			
		||||
@ -114,7 +175,9 @@ constexpr std::tuple<double, double, double, double> marching_squares_lookup(int
 | 
			
		||||
        case 14: return {0, 0.5, 0.5, 0};
 | 
			
		||||
    };
 | 
			
		||||
}
 | 
			
		||||
void write_output_file(std::vector<std::vector<Point>> 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<std::vector<Point>> all_points, const char *f
 | 
			
		||||
    int height = heightmap->height -1;
 | 
			
		||||
    int size = width * height;
 | 
			
		||||
 | 
			
		||||
    for (int j = 0; j < all_points.size(); j++)
 | 
			
		||||
    {
 | 
			
		||||
            OGRFeature *feature;
 | 
			
		||||
            OGRLineString *geometry = new OGRLineString();
 | 
			
		||||
            feature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() );
 | 
			
		||||
            feature->SetField( "Name", j );
 | 
			
		||||
    // gets the lines from the heightmap. Runs the actual algorithm
 | 
			
		||||
    auto lines = create_contours(poLayer->GetLayerDefn(), heightmap, 5);
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
        std::vector<Point> points = all_points[j];
 | 
			
		||||
 | 
			
		||||
        for (int i = 0; i < points.size(); i++) 
 | 
			
		||||
    for (int j = 0; j < lines.size(); j++)
 | 
			
		||||
    {
 | 
			
		||||
 | 
			
		||||
            //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);
 | 
			
		||||
}
 | 
			
		||||
		Loading…
	
	
			
			x
			
			
		
	
		Reference in New Issue
	
	Block a user