2013-11-24 44 views
2

我打算在python中編寫此代碼,但這更像是一個通用算法設計問題,而不是一個與任何特定語言相關的問題。我在寫代碼來模擬一種混合火箭發動機,並且長話短說,這個問題涉及尋找兩者的周長和重疊的圓的很多(可能是成千上萬)構成的圖形的面積。由重疊圓圈構成的區域和周長

我看到了這對stackexchange:

Combined area of overlapping circles

除了我需要找到周邊也。該線程中的某人提到了蒙特卡洛(隨機點猜測)方法,但可以用它來查找周邊以及區域?

在此先感謝。

+0

的區域隨機抽樣方法效果很好,因爲你只關心點的進出身材比例。同樣的事情對於外圍來說也不起作用,因爲太少的點將完全落在它上面,並且「外圍區域」沒有意義。我想你必須根據重疊點來計算周長弧的長度。 – jonrsharpe

+0

http://math.stackexchange.com是您的網站! –

+0

你見過[這](http://stackoverflow.com/questions/13517211/python-how-calculate-a-polygon-perimeter-using-an-osgeo-ogr-geometry-object)的問題?它可能會讓你指向正確的方向。 –

回答

1

比方說,你有圓C1,C2,...,CN。

採取c1和相互交叉的圓圈它。爲c1初始化一個空的圓形扇區列表。

如果交點爲零點或一點,並且c1在另一個圓圈外,則將360度扇區添加到列表中。

如果交點爲零點或一點,並且c1在另一個圓內,則c1的有效輪廓爲0,停止輪廓爲0的處理c1並取下一個圓。

如果交點爲兩點,則將不在另一個圓圈中的c1的扇區添加到c1的圓扇區列表中。

c1的有效提綱是c1扇區列表中所有扇區的交集長度。這個扇區的交集可以由各種不相交的扇區組成。

沖洗和重複的各界人士,然後總結得出的值。

+0

我正在開發一個實現,但是我的數學有點生疏,所以可能需要一段時間。 – Hyperboreus

+0

完成實施,請參閱我的其他答案。 – Hyperboreus

1

我添加了第二個答案,而不是擴展第一個答案,因爲我非常肯定其他答案是正確的,但我不太確定這個答案是否正確。我已經做了一些簡單的測試,似乎可行,但請指出我的錯誤。

這基本上是一個快速的正骯髒執行什麼我前面所說:

from math import atan2, pi 

def intersectLine (a, b): 
    s0, e0, s1, e1 = a [0], a [1], b [0], b [1] 
    s = max (s0, s1) 
    e = min (e0, e1) 
    if e <= s: return (0, 0) 
    return (s, e) 

class SectorSet: 
    def __init__ (self, sectors): 
     self.sectors = sectors [:] 

    def __repr__ (self): 
     return repr (self.sectors) 

    def __iadd__ (self, other): 
     acc = [] 
     for mine in self.sectors: 
      for others in other.sectors: 
       acc.append (intersectLine (mine, others)) 
     self.sectors = [x for x in acc if x [0] != x [1] ] 
     return self 

class Circle: 
    CONTAINS = 0 
    CONTAINED = 1 
    DISJOINT = 2 
    INTERSECT = 3 

    def __init__ (self, x, y, r): 
     self.x = float (x) 
     self.y = float (y) 
     self.r = float (r) 

    def intersect (self, other): 
     a, b, c, d, r0, r1 = self.x, self.y, other.x, other.y, self.r, other.r 
     r0sq, r1sq = r0 ** 2, r1 ** 2 
     Dsq = (c - a) ** 2 + (d - b) ** 2 
     D = Dsq ** .5 
     if D >= r0 + r1: 
      return Circle.DISJOINT, None, None 
     if D <= abs (r0 - r1): 
      return Circle.CONTAINED if r0 < r1 else Circle.CONTAINS, None, None 
     dd = .25 * ((D + r0 + r1) * (D + r0 - r1) * (D - r0 + r1) * (-D + r0 + r1)) ** .5 
     x = (a + c)/2. + (c - a) * (r0sq - r1sq)/2./Dsq 
     x1 = x + 2 * (b - d)/Dsq * dd 
     x2 = x - 2 * (b - d)/Dsq * dd 
     y = (b + d)/2. + (d - b) * (r0sq - r1sq)/2./Dsq 
     y1 = y - 2 * (a - c)/Dsq * dd 
     y2 = y + 2 * (a - c)/Dsq * dd 
     return Circle.INTERSECT, (x1, y1), (x2, y2) 

    def angle (self, point): 
     x0, y0, x, y = self.x, self.y, point [0], point [1] 
     a = atan2 (y - y0, x - x0) + 1.5 * pi 
     if a >= 2 * pi: a -= 2 * pi 
     return a/pi * 180 

    def sector (self, other): 
     typ, i1, i2 = self.intersect (other) 
     if typ == Circle.DISJOINT: return SectorSet ([ (0, 360) ]) 
     if typ == Circle.CONTAINS: return SectorSet ([ (0, 360) ]) 
     if typ == Circle.CONTAINED: return SectorSet ([]) 
     a1 = self.angle (i1) 
     a2 = self.angle (i2) 
     if a1 > a2: return SectorSet ([ (0, a2), (a1, 360) ]) 
     return SectorSet ([ (a1, a2) ]) 

    def outline (self, others): 
     sectors = SectorSet ([ (0, 360) ]) 
     for other in others: 
      sectors += self.sector (other) 
     u = 2 * self.r * pi 
     total = 0 
     for start, end in sectors.sectors: 
      total += (end - start)/360. * u 
     return total 

def outline (circles): 
    total = 0 
    for circle in circles: 
     others = [other for other in circles if circle != other] 
     total += circle.outline (others) 
    return total 

a = Circle (0, 0, 2) 
b = Circle (-2, -1, 1) 
c = Circle (2, -1, 1) 
d = Circle (0, 2, 1) 
print (outline ([a, b, c, d])) 

公式計算從here無恥地竊取兩個圓的交點。