从给定区域采样

Sampling from a given region

我正在探索一种从二维坐标(蓝色区域)上的以下区域进行采样的方法:

当然,基准是我可以抽取随机数(x,y),然后检查它是否与较小的框重叠或在较大的框外。但是,经过一些快速试验,这只会浪费太多的计算资源。

如有任何建议或建议,我们将不胜感激。

可能有一些限制可以允许更简单的解决方案。

所以以下可能不是您的案例的满意解决方案!

但这是一个非常通用的解决方案,这就是为什么我希望在这里 post 没问题。

首先,从绘图上看,矩形似乎总是以原点为中心(两者都是)。如果这是一个有效的假设,则可以简化以下解决方案的部分。

那么,如何使用您建议的 "baseline" 解决方案并不完全清楚。您建议生成点 (x,y),并且对于每个点,您检查它是否包含在内部矩形中。如果它包含在内部矩形中,则将其丢弃。

现在假设您想从蓝色区域采样 100 个点。您必须生成多少点才能确保找到 100 个 被丢弃的点?

这无法确定地解决。或者更正式地说:您不能提供 totally correct 的实现。随机数生成器 可能 总是生成位于内部矩形中的点,因此被丢弃。当然,实际上它不会这样做,但你无法证明这一点,这就是重点。

如果内部矩形与外部矩形相比 "large",这将具有实际意义。您可能需要生成几百万个点才能获得位于内部和外部矩形之间狭窄边缘的 100 个点。


然而,以下是不存在上述问题的解决方案。这是有代价的:它不是一个特别有效的解决方案(尽管如上所述,与基线解决方案相比 "relative efficiency" 取决于矩形的大小和使用模式)。

假设角点的坐标如下图所示:

(x0,y0)                        (x3,y3)
   O------------------------------O
   |                              |
   |  (ix0,iy0)        (ix3,iy3)  |
   |      O----------------O      |
   |      |                |      |
   |      |                |      |
   |      |                |      |
   |      |                |      |
   |      |                |      |
   |      O----------------O      |
   |  (ix1,iy1)        (ix2,iy2)  |
   |                              |
   O------------------------------O
(x1,y1)                        (x2,y2)

(注意坐标是任意的,矩形不一定以原点为中心)

据此,您可以计算可能包含点的区域:

   O------O----------------O------O
   |      |                |      |
   |  R0  |       R1       |  R2  |
   O------O----------------O------|
   |      |                |      |
   |      |                |      |
   |  R2  |                |  R4  |
   |      |                |      |
   |      |                |      |
   O------O----------------O------O
   |  R5  |       R6       |  R7  |
   |      |                |      |
   O------O----------------O------O

现在,当您想对 n 个点进行采样时,您可以为每个点随机 select 这些区域之一,并将该点放置在该区域内的随机位置。

然后需要注意的是该区域的 select 离子:该区域的选择概率应对应于该区域的 面积 ,相对于所有区域的总面积。务实地说,您可以计算所有区域的总面积(如 outer.w*outer.h-inner.w*inner.h),然后计算一个点在区域之一中结束的累积概率分布 R0...R7。从这些累积分布中,您可以将 0.01.0 之间的随机值映射到适当的区域。

这种方法的优点是

  • 它适用于任何矩形(只要外部矩形包含内部矩形...)
  • 你可以直接指定应该生成的点数,预先,它会准确地生成正确的点数
  • 它可以确定性地实现(即完全正确)

下面是显示结果的示例,拖动滑块生成 1...2000 点:

它是使用以下 MCVE 生成的。它是在 Java 中实现的,但是只要你有像 PointRectangle 这样的结构,将它移植到其他语言应该是相当简单的:

import java.awt.BorderLayout;
import java.awt.Color;
import java.awt.Graphics;
import java.awt.Graphics2D;
import java.awt.RenderingHints;
import java.awt.geom.Ellipse2D;
import java.awt.geom.Point2D;
import java.awt.geom.Rectangle2D;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.Random;

import javax.swing.JFrame;
import javax.swing.JPanel;
import javax.swing.JSlider;
import javax.swing.SwingUtilities;

public class RegionNoiseTest
{
    public static void main(String[] args)
    {
        SwingUtilities.invokeLater(() -> createAndShowGUI());
    }

    private static void createAndShowGUI()
    {
        JFrame f = new JFrame();
        f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        f.getContentPane().setLayout(new BorderLayout());

        RegionNoiseTestPanel panel = 
            new RegionNoiseTestPanel();
        f.getContentPane().add(panel, BorderLayout.CENTER);

        JSlider nSlider = new JSlider(1, 2000, 1);
        nSlider.addChangeListener(e -> 
        {
            panel.generatePoints(nSlider.getValue());
        });
        nSlider.setValue(100);
        f.getContentPane().add(nSlider, BorderLayout.SOUTH);

        f.setSize(500,450);
        f.setLocationRelativeTo(null);
        f.setVisible(true);
    }

}

class RegionNoiseTestPanel extends JPanel
{
    private final Rectangle2D outer;
    private final Rectangle2D inner;
    private List<Point2D> points;


    RegionNoiseTestPanel()
    {
        outer = new Rectangle2D.Double(50, 50, 400, 300);
        inner = new Rectangle2D.Double(90, 100, 300, 200);
    }

    public void generatePoints(int n)
    {
        this.points = createPoints(n, outer, inner, new Random(0));
        repaint();
    }

    @Override
    protected void paintComponent(Graphics gr)
    {
        super.paintComponent(gr);
        Graphics2D g = (Graphics2D)gr;
        g.setRenderingHint(
            RenderingHints.KEY_ANTIALIASING, 
            RenderingHints.VALUE_ANTIALIAS_ON);

        g.setColor(new Color(220, 220, 220));
        g.fill(outer);
        g.setColor(new Color(160, 160, 160));
        g.fill(inner);

        if (points != null)
        {
            g.setColor(Color.BLUE);
            for (Point2D p : points)
            {
                double r = 2;
                double x = p.getX();
                double y = p.getY();
                g.fill(new Ellipse2D.Double(x - r, y - r, r + r, r + r));
            }
        }


    }

    private static List<Point2D> createPoints(
        int n, Rectangle2D outer, Rectangle2D inner, Random random)
    {
        List<Rectangle2D> regions = computeRegions(outer, inner);
        double cumulativeRegionAreas[] = new double[8];
        double outerArea = outer.getWidth() * outer.getHeight();
        double innerArea = inner.getWidth() * inner.getHeight();
        double relevantArea = outerArea - innerArea;
        double areaSum = 0;
        for (int i = 0; i < regions.size(); i++)
        {
            Rectangle2D region = regions.get(i);
            double area = region.getWidth() * region.getHeight();
            areaSum += area;
            cumulativeRegionAreas[i] = areaSum / relevantArea;
        }

        List<Point2D> points = new ArrayList<Point2D>();
        for (int i=0; i<n; i++)
        {
            points.add(createPoint(
                regions, cumulativeRegionAreas, random));
        }
        return points;
    }

    private static List<Rectangle2D> computeRegions(
        Rectangle2D outer, Rectangle2D inner)
    {
        List<Rectangle2D> regions = new ArrayList<Rectangle2D>();
        for (int r = 0; r < 8; r++)
        {
            regions.add(createRegion(outer, inner, r));
        }
        return regions;
    }

    private static Point2D createPoint(
        List<Rectangle2D> regions, 
        double normalizedCumulativeRegionAreas[], 
        Random random)
    {
        double alpha = random.nextDouble();
        int index = Arrays.binarySearch(normalizedCumulativeRegionAreas, alpha);
        if (index < 0)
        {
            index = -(index + 1);
        }
        Rectangle2D region = regions.get(index);
        double minX = region.getMinX();
        double minY = region.getMinY();
        double maxX = region.getMaxX();
        double maxY = region.getMaxY();
        double x = minX + random.nextDouble() * (maxX - minX);
        double y = minY + random.nextDouble() * (maxY - minY);
        return new Point2D.Double(x, y);
    }

    private static Rectangle2D createRegion(
        Rectangle2D outer, Rectangle2D inner, int region)
    {
        double minX = 0;
        double minY = 0;
        double maxX = 0;
        double maxY = 0;
        switch (region) 
        {
            case 0:
                minX = outer.getMinX();
                minY = outer.getMinY();
                maxX = inner.getMinX();
                maxY = inner.getMinY();
                break;

            case 1:
                minX = inner.getMinX();
                minY = outer.getMinY();
                maxX = inner.getMaxX();
                maxY = inner.getMinY();
                break;

            case 2:
                minX = inner.getMaxX();
                minY = outer.getMinY();
                maxX = outer.getMaxX();
                maxY = inner.getMinY();
                break;

            case 3:
                minX = outer.getMinX();
                minY = inner.getMinY();
                maxX = inner.getMinX();
                maxY = inner.getMaxY();
                break;

            case 4:
                minX = inner.getMaxX();
                minY = inner.getMinY();
                maxX = outer.getMaxX();
                maxY = inner.getMaxY();
                break;

            case 5:
                minX = outer.getMinX();
                minY = inner.getMaxY();
                maxX = inner.getMinX();
                maxY = outer.getMaxY();
                break;

            case 6:
                minX = inner.getMinX();
                minY = inner.getMaxY();
                maxX = inner.getMaxX();
                maxY = outer.getMaxY();
                break;

            case 7:
                minX = inner.getMaxX();
                minY = inner.getMaxY();
                maxX = outer.getMaxX();
                maxY = outer.getMaxY();
                break;
        }
        return new Rectangle2D.Double(minX, minY, maxX - minX, maxY - minY);
    }

}

我仍然很好奇是否有人找到了一种优雅的、确定性的方法 而不是 需要为要生成的点定义各种 "regions"。 ..

如果蓝色区域与原点对称,您可以从从单位正方形采样的随机点创建映射。考虑以下伪代码函数,假设两个矩形都以原点为中心:

def sample():
    sample point x_base from [-1, 1] and calculate x = sign(x_base)*x1 + x_base*(x2-x1)
    sample point y_base from [-1, 1] and calculate y = sign(y_base)*y1 + y_base*(y2-y1)
    if (x,y) == (0,0) then:
        # recursively sample again in the rare case of sampling 0 for both dimensions
        return sample()
    else:
        return (x,y)

编辑: 正如 Marco13 所指出的那样,此解决方案无法从整个蓝色区域正确采样。查看他的回答以获得更好的方法。

这里的解决方案假设蓝色区域是对称的并且以原点为中心,所以有4个参数(x1, x2, y1, y2)。想象一下,在蓝色区域内是另一个具有相同比例但按比例缩小的区域,以便该其他区域的外边界恰好适合蓝色区域的内边界。如果我们随机生成一个点并且它位于这个内部区域内,我们可以通过将 x 和 y 分别缩放 x2/x1 和 y2/y1 将其映射到蓝色区域中的一个点。现在想象一下这个区域内部的另一个区域,以及它内部的另一个区域,无穷无尽。任何点(原点除外)都可以通过将其放大适当的次数来映射到蓝色区域中的一个点:

// generate a random point:
double x = 0.0, y = 0.0;
while(x == 0.0 && y == 0.0) // exclude the origin
{
    x = random.NextDouble() * x2;
    y = random.NextDouble() * y2;
}

// map to the blue region
while(x < x1 && y < y1)
{
    x *= (x2 / x1);
    y *= (y2 / y1);
}

// randomly choose a quadrant:
int r = random.Next(0, 4);
if((r & 1) != 0)
    x = -x;
if((r & 2) != 0)
    y = -y;

但是由于第二个 while 循环(实际上保证第一个 while 循环永远不会 运行 超过一次),这并不是很好。可以使用对数消除循环:

// map to the blue region
if(x < x1 && y < y1)
{
    double xpower = Math.Ceiling((Math.Log(x1) - Math.Log(x)) / Math.Log(x2/x1));
    double ypower = Math.Ceiling((Math.Log(y1) - Math.Log(y)) / Math.Log(y2/y1));
    double power = Math.Min(xpower, ypower);
    x *= Math.Pow(x2/x1, power);
    y *= Math.Pow(y2/y1, power);
}