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

代碼在每一章節最後

 

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

咱們知道生成矩形內的隨機點比較容易，只要分別隨機生成相應的橫座標和縱座標，好比隨機生成範圍[-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()

---

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

算法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