如何在python中栅格化一批多边形

Question

我有一批多边形，它们可以是 NumPy 数组、火炬张量或任何其他形状为 [230_000,3,2] 的 nd 数组。 230_000 是多边形的个数，3 表示它是一个三角形，2 是它的 x 和 y 坐标。

我还有一批形状为[230_000,3,3即颜色]的特征，最好在光栅化时插值特征的颜色。

我正在尝试找到一种方法来栅格化所有这些，以使输出看起来像这样：

算法的规则不多。最重要的是在 nd 阵列上栅格化的填充多边形并且速度很快。批量大小将是数十万个多边形。 for 循环不会有帮助。

Answer 1

如果您打算使用 python 创建栅格，我会查看 Rasterio。如果您的主要目标是使用 Python 创建上述可视化，那么更简单的方法是使用 shapely 和 matplotlib。看起来你可能会丢失一块——你还需要一个生成形状的函数。你提到：

230_000 是多边形的个数，3 表示它是一个三角形，2 是它的 x 和 y 坐标

2 是什么点的x 和y 坐标？了解您需要生成的形状类型需要更多信息。

AFAIK，你需要做的是：

• 迭代你的数组（没办法）
• 定义一个函数，该函数接受一个单独的数组并从提供的输入https://autogis-site.readthedocs.io/en/latest/notebooks/L1/geometric-objects.html 返回一个匀称的多边形 从那里我会：
• 转换为geopandas数据框https://geopandas.org/en/stable/docs/reference/geodataframe.html
• 将 GeoDataFrame 转换为 GeoJSON https://geopandas.org/en/stable/docs/reference/api/geopandas.GeoDataFrame.to_json.html
• 将geojson传递给rasterio https://rasterio.readthedocs.io/en/latest/cli.html

Answer 2

我想出了一种方法，可以使用仅张量变换来栅格化一批三角形。我的实现在 PyTorch 中，但在其他库中会类似。它不是超级有效，但效果很好。该算法动态计算具有不同长度的像素和索引，从而减少较小三角形的计算量。源码如下：

import torch
import torch.nn as nn
import cv2

class TriRasterizer(nn.Module):
    def __init__(self,shape=(480,640),channels:int=3,_trainable=False):
        super(TriRasterizer, self).__init__()
        assert len(shape) == 2
        self._shape = shape
        self._backdrop = nn.Parameter(torch.zeros(*shape,channels,dtype=torch.float),requires_grad=_trainable)
        # Feel free to mess around with the quality
        self.quality = 0.1
        ...

    @staticmethod
    @torch.jit.script
    def _fill(poly,area,_type:torch.dtype=torch.float32):
        """
        The main part
        The algorithm takes the tensor points, repeats interleave them by their area
        Then for each tiangles in the batch:
            there should be the-area more triangles exactly the same
            the goal is for each of those extra triangles,
                interpolate each of those points such that: each points fills up the traingle

        """
        # calculate center of the base
        base_p = (poly[:,0] + poly[:,1])/2
        a_diff = (poly[:,2]-base_p) . detach()

        # then height
        heights = torch.sqrt_(torch.square_(a_diff[...,0]) + torch.square_(a_diff[...,1]))

        # comulative sum of the areas
        # This gives positional embeding for each points when subtracted from torch.arange
        cum   = torch.cumsum(area,dim=0).int()
        poly  = poly.type(_type)
        heights  = heights.type(torch.int16)
        poly    = poly.repeat_interleave(area,dim=0)
        cum_    = cum.repeat_interleave(area,dim=0)
        heights = heights.repeat_interleave(area,dim=0)

        # normalized makes the first point = 0 and last = 1 for the repeated points
        normalized = 1-(cum_-torch.arange(len(cum_)))/cum_

        # this is the exponential, it reduces pixels as it gets closer to the tip of the triangle
        exp = 1-torch.sqrt_(1-normalized)
        exp = exp.type(torch.float32)

        # This iterpolates more of the x as it gets to the tip of the triangle
        ij__interp = (exp % (torch.tensor(1,dtype=torch.float16)/heights))[...,None]
        # * heights[...,None] just brings the exponentials max back to 1
        # the random normal washes the points and makes it appear blurry
        ij_interp = ij__interp*heights[...,None] #+ torch.randn_like(ij__interp)/70

        # this makes the makes the c interpolation exponentialy faster
        k_interp = (exp[...,None] - ij__interp)  #+ torch.randn_like(ij__interp)/70

        # Finally interpolate them
        # weigh a and b first with respect to the calculated exponential
        ij_trunc = poly[:,0] * (1-ij_interp) + poly[:,1]*ij_interp
        # then weigh the output with c with respect to the modulus of the exponential
        filled   = poly[:,2] * k_interp  + ij_trunc * (1-k_interp)
        return filled

    def _compute_area(self,poly):
        # Formula to find area from a set of points
        area = torch.abs(0.5 * ((poly[..., 0, 0] * (poly[..., 1, 1] - poly[..., 2, 1])) +
                   (poly[...,1,0]*(poly[...,2,1]-poly[...,0,1])) +
                   (poly[...,2,0]*(poly[...,0,1]-poly[...,1,1]))))
        return area

    def forward(self,traingles,features):
        return self.rasterize(traingles,features)


    @torch.no_grad()
    def rasterize(self,traingles,features):
        raster = self._backdrop.clone()
        # first find the area
        # the lower the quality the more spaces in between pixels
        area = (self._compute_area(traingles)*self.quality).int()

        # fill in points based on computed area
        filled_cords = self._fill(traingles,area,_type=torch.int16)
        filled_cords = filled_cords.long()

        # clamp the points to prevent out of bounds exception
        filled_cords[...,0] = filled_cords[...,0].clamp(min=0,max=self._shape[1]-1)
        filled_cords[...,1] = filled_cords[...,1].clamp(min=0,max=self._shape[0]-1)

        # also fill in features based on computed area
        filled_feats = self._fill(features,area,_type=torch.float32)
        pixels = (filled_cords[...,1],filled_cords[...,0])

        # compute the times the same point has been placed on the raster
        sum_ = torch.zeros_like(raster).index_put(pixels,torch.ones_like(filled_feats),accumulate=True)

        # rasterize them using indexing
        raster = raster.index_put(pixels,filled_feats,accumulate=False)

        # average those sums
        raster = raster/sum_
        return raster
#

if __name__ == '__main__':
    tiangles = torch.tensor([
        [[100, 0], [0, 200], [20, 300]]
    ], dtype=torch.float) + 100

    print(tiangles.shape)
    tiangles = tiangles.expand(3,-1,-1).clone()

    features = torch.randn(*tiangles.shape[:-1], 3)
    rasterizer = TriRasterizer()
    j = torch.zeros(*tiangles.shape)

    while True:
        features = features*0.01 + torch.randn(*tiangles.shape[:-1], 3)*0.09
        j += torch.randn(*tiangles.shape)*0.1
        tiangles += j
        raster = rasterizer(tiangles,features+1)
        cv2.imshow('x',raster.numpy())
        cv2.waitKey(20)