How can we identify and delineate holes in a 2D point set representing soil sample locations?

Finding Holes in 2D Point Sets

The task is to find the holes in a set of 2D points within a cartesian grid system. The points represent soil sample locations, and holes could include giant rocks, swampy places, or lakes/ponds. The goal is to find the concave polygon that roughly defines these areas, adjusting the sensitivity of the algorithm to control the roughness or smoothness of the polygon.

Solution Approach

Steps:

1. Create a density map: Convert the point set to a bitmap or 2D array by scaling and projecting each point onto a grid. Calculate the density (number of points) for each cell.
2. Identify holes: Find cells with zero density or below a given threshold.
3. Segment hole areas: Create horizontal and vertical lines covering these holes, grouping them by proximity to form hole segments.
4. Polygonize hole segments: Convert the segments into concave polygons. Sort the points to ensure proper connectivity and remove duplicates.

Example Implementation (C#):

using System;
using System.Collections.Generic;

public class Holes
{
    // Density map (2D array)
    private int[][] map;

    // List of hole segments (lines)
    private List<Line> segments;

    // Polygonized holes (concave polygons)
    private List<Polygon> holes;

    // Polygonization tolerance (higher value = smoother polygons)
    private double tolerance;

    // Initializes the hole detection algorithm.
    public Holes(int[][] points, int mapSize, double tolerance)
    {
        if (points == null || mapSize <= 0 || tolerance <= 0)
        {
            throw new ArgumentException("Invalid arguments");
        }
        
        // Initialize the variables
        this.map = new int[mapSize][mapSize];
        this.tolerance = tolerance;
        this.segments = new List<Line>();
        this.holes = new List<Polygon>();
        
        // Create density map
        CreateDensityMap(points, mapSize);
    }

    // Identifies holes in the density map.
    public void FindHoles()
    {
        if (map == null || map.Length == 0)
        {
            throw new InvalidOperationException("Density map not initialized.");
        }
        
        // Find hole cells
        List<Cell> holeCells = FindCells(0);
        
        // Group hole cells into segments
        List<List<Line>> lineGroups = GroupLines(holeCells);
        
        // Polygonize segments
        PolygonizeSegments(lineGroups);
    }

    // Helper functions for hole detection.

    private void CreateDensityMap(int[][] points, int mapSize)
    {
        // Scale and project points onto a grid
        for (int i = 0; i < points.Length; i++)
        {
            double scaledX = points[i][0] / points[0][0] * mapSize;
            double scaledY = points[i][1] / points[0][1] * mapSize;
            int x = (int)scaledX;
            int y = (int)scaledY;
            
            // Increment count in density map
            map[x][y]++;
        }
    }

    private List<Cell> FindCells(int threshold)
    {
        List<Cell> holeCells = new List<Cell>();
        
        for (int i = 0; i < map.Length; i++)
        {
            for (int j = 0; j < map[i].Length; j++)
            {
                if (map[i][j] == 0 || map[i][j] <= threshold)
                {
                    holeCells.Add(new Cell(i, j));
                }
            }
        }
        
        return holeCells;
    }

    private List<List<Line>> GroupLines(List<Cell> holeCells)
    {
        // Group lines by proximity
        List<List<Line>> lineGroups = new List<List<Line>>();
        foreach (Cell holeCell in holeCells)
        {
            List<Line> group = null;
            
            // Find existing group or create a new one
            for (int i = 0; i < lineGroups.Count; i++)
            {
                if (lineGroups[i].Find(line => line.Proximity(holeCell) <= tolerance) != null)
                {
                    group = lineGroups[i];
                    break;
                }
            }
            
            if (group == null)
            {
                group = new List<Line>();
                lineGroups.Add(group);
            }
            
            // Add horizontal/vertical lines
            group.Add(new Line(holeCell.x, holeCell.y, true));
            group.Add(new Line(holeCell.x, holeCell.y, false));
        }
        
        return lineGroups;
    }

    private void PolygonizeSegments(List<List<Line>> lineGroups)
    {
        foreach (List<Line> lineGroup in lineGroups)
        {
            Polygon polygon = PolygonizeSegment(lineGroup);
            if (polygon != null)
            {
                holes.Add(polygon);
            }
        }
    }

    private Polygon PolygonizeSegment(List<Line> lineSegment)
    {
        // Sort lines by angle (convex hull algorithm)
        lineSegment.Sort((a, b) => a.Angle.CompareTo(b.Angle));
        
        // Remove duplicate lines
        List<Line> uniqueLines = new List<Line>();
        foreach (Line line in lineSegment)
        {
            if (uniqueLines.Count == 0 || uniqueLines[uniqueLines.Count - 1].Angle != line.Angle)
            {
                uniqueLines.Add(line);
            }
        }
        
        // Polygonize lines
        List<Point> points = new List<Point>();
        for (int i = 0; i < uniqueLines.Count; i++)
        {
            Point point = null;
            Line currentLine = uniqueLines[i];
            
            if (uniqueLines[(i + 1) % uniqueLines.Count].Angle - currentLine.Angle > Math.PI)
            {
                point = currentLine.GetIntersection(uniqueLines[(i + 1) % uniqueLines.Count], true);
            }
            else
            {
                point = currentLine.GetIntersection(uniqueLines[(i + 1) % uniqueLines.Count], false);
            }
            
            if (point != null)
            {
                points.Add(point);
            }
        }
        
        return new Polygon(points);
    }

    // Helper classes for line/polygon representation.

    private class Line
    {
        public int x1, y1, x2, y2;
        public double angle;
        public bool isHorizontal;

        public Line(int x, int y, bool isHorizontal)
        {
            if (isHorizontal)
            {
                x1 = 0; y1 = y;
                x2 = map.GetLength(0) - 1; y2 = y;
            }
            else
            {
                x1 = x; y1 = 0;
                x2 = x; y2 = map[0].GetLength(0) - 1;
            }
            
            this.angle = Math.Atan2(y2 - y1, x2 - x1);
            this.isHorizontal = isHorizontal;
        }

        public double Angle { get { return angle; } }

        public double Proximity(Cell cell)
        {
            double distX, distY;
            if (isHorizontal)
            {
                distX = cell.x - x1;
                distY = cell.y - y1;
            }
            else
            {
                distX = cell.x - x2;
                distY = cell.y - y2;
            }
            
            return Math.Sqrt(distX * distX + distY * distY);
        }

        public Point GetIntersection(Line other, bool isConvex)
        {
            double denominator, numerator, tx, ty;
            
            if (isHorizontal)
            {
                denominator = (other.y2 - other.y1) - (y2 - y1);
                numerator = ((other.x2 - other.x1) * (y1 - other.y1)) - ((x2 - x1) * (other.y2 - other.y1));
                tx = numerator / denominator;
                ty = other.y1 + ((tx - other.x1) * (other.y2 - other.y1)) / (other.x2 - other.x1);
            }
            else
            {
                denominator = (other.x2 - other.x1) - (x2 - x1);

The above is the detailed content of How can we identify and delineate holes in a 2D point set representing soil sample locations?. For more information, please follow other related articles on the PHP Chinese website!