均勻的生成圓和三角形內的隨機點

代碼在每一章節最后

 

一、均勻生成圓內的隨機點

我們知道生成矩形內的隨機點比較容易，只要分別隨機生成相應的橫坐標和縱坐標，比如隨機生成范圍[-10,10]內橫坐標x，隨機生成范圍[-20,20]內的縱坐標y，那么(x,y)就是生成的隨機點。由此，我們很容易的想到了算法1

算法1（正確的）：

每個圓對應一個外切矩形，我們隨機生成矩形內的點，如果該點在圓內，就返回改點，否則重新生成直到生成的點在圓內。

該方法的缺點是有可能連續幾次都生成不了符合要求的點。（可以求得：每次生成點時，該點有 的概率在圓內）

 

算法2（錯誤的）：

可能有的人想到該方法，根據圓的解析式 （假設圓心在原點）我們可以先隨機生成[-R, R]范圍內橫坐標x，然后生成 范圍內的隨機數y，(x,y)就是需要的點。

我們寫程序模擬了該過程，從下圖可以看出，我們可以看到當x靠近圓的邊緣使，y的范圍減小，因此兩邊邊緣的點較密集，靠近圓心的點較稀疏。

 

算法3（錯誤的）：

還可以根據極坐標，圓內的點可以如下描述

x = r*sin(theta)

y = r*cos(theta)

其中0 <= r <= R, 0 <= theta < 360

先隨機生成[0, 360)內的theta，然后隨機生成[0, R]內的r。

theta固定后，當r越靠近R，即點越靠近圓的邊緣，因此如果r是[0，R]等概率生成的，那么圓的邊緣的點會比靠近圓心處要稀疏。

算法4（正確的）：

和算法3一樣還是根據極坐標

x = r*sin(theta)

y = r*cos(theta)

其中0 <= r <= R, 0 <= theta < 360                      本文地址

先隨機生成[0, 360)內的theta，然后隨機生成[0, 1]內的k, 且r = sqrt(k)*R。

根據根號函數的性質，這樣使得r的分布在k靠近1（靠近邊緣）的地方點變得較密，具體證明可參考here, 也可以參考論文“矩形和橢圓內均勻分布隨機點定理及應用”，圓是橢圓的特例

以上的4個算法的代碼如下(Python3.3):

import numpy as np
import matplotlib.pyplot as plt
import random
import math

#博客算法1
def GeneratePointInCycle1(point_num, radius):
    for i in range(1, point_num+1):
        while True:
            x = random.uniform(-radius, radius)
            y = random.uniform(-radius, radius)
            if (x**2) + (y**2) < (radius**2):
                break
        plt.plot(x, y, '*', color = "black")

#博客算法2
def GeneratePointInCycle2(point_num, radius):
    for i in range(1, point_num+1):
        x = random.uniform(-radius, radius)
        y_max = math.sqrt(radius**2 - x**2)
        y = random.uniform(-y_max, y_max)
        plt.plot(x, y, '*', color = "black")

#博客算法3
def GeneratePointInCycle3(point_num, radius):
    for i in range(1, point_num+1):
        theta = random.random()*2*pi;
        r = random.uniform(0, radius)
        x = r*math.sin(theta)
        y = r*math.cos(theta)
        plt.plot(x, y, '*', color = "black")

#博客算法4
def GeneratePointInCycle4(point_num, radius):
    for i in range(1, point_num+1):
        theta = random.random()*2*pi;
        r = random.uniform(0, radius)
        x = math.sin(theta)* (r**0.5)
        y = math.cos(theta)* (r**0.5)
        plt.plot(x, y, '*', color = "black")      


pi = np.pi
theta = np.linspace(0, pi*2, 1000)
R = 1
x = np.sin(theta)*R
y = np.cos(theta)*R

plt.figure(figsize=(6,6))
plt.plot(x,y,label = "cycle",color="red",linewidth=2)
plt.title("cycyle")
GeneratePointInCycle4(4000, R) #修改此處來顯示不同算法的效果
plt.legend()
plt.show()

---

一、均勻生成三角形內的隨機點

算法5（錯誤的）

對於三角形ABC和一點P，可以有如下的向量表示：

p點在三角形內部的充分必要條件是：1 >= u >= 0, 1 >= v >= 0, u+v <= 1。

先生成[0,1]的隨機數u，然后生成[0, 1-u]內的隨機數v，u、v生成后，就可以得到p點的坐標:

由下圖可知，該算法生成的點在靠近A點處較濃密

 

算法6（正確的）

如圖所示，三角形ABC有與之對應的矩形ABNM，且矩形面積是三角形的兩倍，三角形ADC和CMA全等，CDB和BNC全等。

我們可以先生成矩形ABNM內的隨機點P，如果P剛好在三角形ABC中，那么符合要求；如果P不在三角形ABC中，P要么在AMC中，要么在BNC中，如圖P在BNC中，我們求P關於BC中點的的中心對稱點，該點一定在三角形中。P在AMC中同理。這樣可以保重三角形外的點都可以均勻的一一對應到三角形內部。

后面的代碼中，為了簡化計算，我們假設AB是平行X軸的。

對於生成任意多邊形內的隨機點，我們可以把它分割成三角形，再來生成隨機點。

算法5和算法6的代碼如下(Python3.3):

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.lines as line
import random
import math

#對應博客算法5
def GeneratePointInTriangle1(point_num, pointA, pointB, pointC):
    for i in range(1, point_num+1):
        u = random.uniform(0.0, 1.0)
        v = random.uniform(0.0, 1.0 - u)
        pointP = u*pointC + v*pointB + (1.0-u-v)*pointA;
        plt.plot(pointP[0], pointP[1], '*', color = "black")

#根據向量叉乘計算三角形面積，參考 http://www.cnblogs.com/TenosDoIt/p/4024413.html
def ComputeTriangleArea(pointA, pointB, pointC):
    return math.fabs(np.cross(pointB - pointA, pointB - pointC)) / 2.0

#判斷點P是否在三角形ABC內,參考 http://www.cnblogs.com/TenosDoIt/p/4024413.html
def IsPointInTriangle(pointA, pointB, pointC, pointP):
    area_abc = ComputeTriangleArea(pointA, pointB, pointC)
    area_pab = ComputeTriangleArea(pointA, pointB, pointP)
    area_pbc = ComputeTriangleArea(pointP, pointB, pointC)
    area_pac = ComputeTriangleArea(pointP, pointA, pointC)
    return math.fabs(area_pab + area_pac + area_pbc - area_abc) < 0.000001

#計算一個點關於某一點的中心對稱點
def ComputeCentralSymmetryPoint(point_src, point_center):
    return np.array([point_center[0]*2-point_src[0], point_center[1]*2-point_src[1]])

#對應博客算法6
def GeneratePointInTriangle2(point_num, pointA, pointB, pointC):
    for i in range(1, point_num+1):
        pointP = np.array([random.uniform(pointA[0], pointB[0]), random.uniform(pointA[1], pointC[1])])
        if not IsPointInTriangle(pointA, pointB, pointC, pointP):
            if pointP[0] > pointC[0]:
                pointP = ComputeCentralSymmetryPoint(pointP, np.array([(pointC[0] + pointB[0])/2, (pointC[1] + pointB[1])/2]))
            else:
                pointP = ComputeCentralSymmetryPoint(pointP, np.array([(pointC[0] + pointA[0])/2, (pointC[1] + pointA[1])/2]))
        plt.plot(pointP[0], pointP[1], '*', color = "black")        


fig = plt.figure()
#三角形三個頂點
pointA = np.array([0,1])
pointB = np.array([3,1])
pointC = np.array([1,2])

plt.plot([pointA[0],pointB[0]], [pointA[1],pointB[1]])
plt.plot([pointA[0],pointC[0]], [pointA[1],pointC[1]])
plt.plot([pointB[0],pointC[0]], [pointB[1],pointC[1]])

GeneratePointInTriangle2(1500, pointA, pointB, pointC) #修改此處來顯示不同算法的效果

plt.ylim(0.5,2)
plt.xlim(0,3)
plt.show()

【版權聲明】轉載請注明出處：http://www.cnblogs.com/TenosDoIt/p/4025221.html