Python复杂形状矩阵处理：多边形区域元素均值计算问题

Hey there! Let’s work through how to solve this problem—this is a common spatial computing task, and I’ve got a straightforward, reliable approach for you that avoids the complexity of trying to transform the polygon into a square.

Core Idea #

To get the average, we just need to:

• Identify every matrix cell that lies inside or on the boundary of your polygon.
• Sum the values of those cells.
• Divide the total sum by the number of qualifying cells.

The key challenge here is accurately checking whether a matrix cell belongs to the polygon region. Let’s break this down step by step.

Step 1: Clarify Coordinate Mapping #

First, make sure you’re clear on how your matrix indices map to the x/y coordinates of your polygon points. For example:

• Does matrix row i, column j correspond to a cell with center at (j + 0.5, i + 0.5)? (This is standard for grid-based systems, where each cell spans from integer coordinates.)
• Double-check that your polygon points use the same coordinate system as your matrix.

Step 2: Implement Point-in-Polygon Check #

The go-to method for this is the ray casting algorithm—it’s efficient and easy to implement. Here’s how it works:
Shoot a horizontal ray (e.g., straight to the right) from the center of your matrix cell. Count how many times this ray crosses the edges of the polygon. If the count is odd, the cell is inside the polygon; if even, it’s outside. We’ll also handle edge cases where the ray passes through a vertex or lies exactly on an edge.

Here’s a sample implementation (adapt this to your language):

def is_inside(point, polygon):
    x, y = point
    inside = False
    num_vertices = len(polygon)
    
    for i in range(num_vertices):
        # Get the current edge's start and end points
        j = (i + 1) % num_vertices
        x1, y1 = polygon[i]
        x2, y2 = polygon[j]
        
        # Check if the point's y-coordinate is within the edge's y-range
        if ((y1 > y) != (y2 > y)):
            # Calculate the x-coordinate where the ray intersects the edge
            x_intersect = ((y - y1) * (x2 - x1)) / (y2 - y1) + x1
            # If the intersection is to the right of the point, flip the inside flag
            if x <= x_intersect:
                inside = not inside
    
    # Optional: Add check for points exactly on edges if needed
    return inside

Step 3: Traverse the Matrix and Compute the Average #

Now loop through every cell in your matrix, check if its center is inside the polygon, and track the sum and count of qualifying cells.

Example code (assuming your matrix is a 2D list called matrix, and your polygon points are stored in polygon_points):

total_sum = 0.0
cell_count = 0
rows = len(matrix)
cols = len(matrix[0]) if rows > 0 else 0

for i in range(rows):
    for j in range(cols):
        # Get the center of the current matrix cell (adjust this if your coordinate system differs)
        cell_center = (j + 0.5, i + 0.5)
        if is_inside(cell_center, polygon_points):
            total_sum += matrix[i][j]
            cell_count += 1

# Handle edge case where no cells are inside the polygon
average = total_sum / cell_count if cell_count > 0 else 0.0

Step 4: Optional Optimization #

If your matrix is very large, you can speed things up by first finding the bounding box of your polygon (the smallest rectangle that contains the polygon). Then only iterate over matrix cells that fall within this bounding box—this reduces the number of point-in-polygon checks you need to run.

Why Your Parameterization Approach Might Have Failed #

Transforming the polygon into a square (using affine transformations, for example) is technically possible, but it introduces a lot of complexity:

• You’d need to ensure the transformation is invertible.
• Discrete matrix cells don’t map cleanly to transformed coordinates, requiring interpolation to estimate values in the square.
• You’d have to map the transformed square back to the original matrix, which adds more room for error.

The ray casting method is simpler, more direct, and avoids all these issues.

内容的提问来源于stack exchange，提问作者Doron bs