如何生成多個互不重疊的不同半徑圓？
如圖所示，採用什麼樣的思路處理才能防止隨機生成的圓之間產生重疊。另補充一個條件，即，圓需按照一定半徑範圍生成，且能控制圓形面積與整體四邊形面積的比例，（如圓半徑範圍為1-2mm,圓形面積與整體面積之比為0.5)。此外，筆者還有一個疑問，如何對圓與圓之間是否產生接觸進行判斷？
這個問題就算沒人邀請，我也一定要答。有一次和公司一個美術扯皮，最後他竟然威脅我說：「我畫個圈圈詛咒你！」我立馬反擊道：「我寫個程序，畫一大坨圈圈詛咒你！而且我還會用N種方法畫。」
（1）第一種方法
寫個循環，循環中隨機生成一個圈圈的圓心坐標位置、半徑。然後判斷圈圈是否與已有的圈圈發生碰撞。如果有碰撞讓它滾蛋，否則畫圈圈。（2）第二種方法
還是寫個循環，循環中隨機生成一個圈圈的圓心坐標位置，然後判斷這個坐標點是否在已有的圈圈內部，如果是就讓它滾蛋，否則計算出與它到距離最近圈圈的距離，並以此畫圈圈。（3）第三種方法
有時需要將一堆指定半徑的圈圈，無碰撞的放置在某一空間中。這就需要一套物理碰撞演算法，感覺很難，不要怕，圈圈之間的碰撞檢測是所有碰撞檢測中最容易的。碰到這一情況，首先要找到一個中心點，然後將所有的圈圈按照與該中心點的距離遠近進行排序，依次將每一個圈圈朝向黨中心靠近，最後實現大和諧。這種方法中需要克服的一個難點是，如何實現三個圈圈的兩兩相碰。解決這一問題，需要用到一個數學公式：餘弦定理，已知三角形的三邊長度，求三個角度。這算是我在實際編程中用過最高深的數學知識之一。我寫的一個小遊戲：連泡泡用的就是這一方法。
畫圈圈，就這麼簡單搞定。
發一下該演算法的C代碼：inline void hash22(float px, float py, float hx, float hy)
{
float n = sinf(px * 7.0f + py * 157.0f);
hx = 2097152.0f*n;
hy = 262144.0f*n;
hx = hx - floorf(hx);
hy = hy - floorf(hy);
}
// 查找最近的相鄰頂點
float nearest_dis_2d(float x, float y)
{
float ox, oy;
hash22(x, y, ox, oy); // 當前頂點偏移
float d, r = 2.0f;
float vx, vy;
for(int i = -1; i &<= 1; i++)
{
for(int j = -1; j &<= 1; j++)
{
if (i == 0 j == 0)
{
continue;
}
hash22(i + x, j + y, vx, vy);
vx += float(i) - ox;
vy += float(j) - oy;
d = vx*vx + vy*vy;
if (d &< r)
{
r = d;
}
}
}
return sqrtf(r);
}
// 使用worley畫圈圈
float circle_worley2d(float x, float y)
{
float fx = floorf(x);
float fy = floorf(y);
float gx = x - fx;
float gy = y - fy;
float vx, vy;
float d, r = 3.0f;
float si, sj;
for(int i = -1; i &<= 1; i++)
{
for(int j = -1; j &<= 1; j++)
{
hash22(i + fx, j + fy, vx, vy);
vx += float(i) - gx;
vy += float(j) - gy;
d = vx*vx + vy*vy;
if (d &< r)
{
r = d;
si = i + fx;
sj = j + fy;
}
}
}
d = sqrtf(r) / nearest_dis_2d(si, sj) * 2.0f;
return d;
}
代碼中調用函數circle_worley2d，輸入是平面上的一個二維坐標，返回值為1時，表示該坐標位於圓上。用上面的代碼可以生成如下圖像：
（6）第六種方法:
上一種方法中,float circle_worley2d(float x, float y)函數,可以很容易地改成三維的float circle_worley3d(float x, float y, float z)。
用三維的可以生成空間中若干個互相不重疊的球,如下圖所示:
當然,擴展到三維後的circle_worley3d函數，可以在任意一個三維圖形上畫圈圈。下圖為將許多圈圈畫在一個圈圈上:
看到這個問題感覺特別感興趣，於是就照著 @葉飛影 和 @柳傾 提到的物理模擬的方法編了個MATLAB模擬程序，做了點微小的貢獻。模擬的結果大致是這樣的：
動圖鏈接：Imgur: The most awesome images on the InternetRemarks:
- 用了很簡單的spring-damper-mass模型。受 @葉飛影 的啟發，所有小球都受到向中心的引力。
- 這個簡單模型帶來的問題就是：所有的球或多或少都會有些形變，而且中心的球的形變尤為嚴重一些。假如要減小形變就得增大彈性係數，但那樣又會出現各種各樣的數值問題……球的半徑大小差距太大也會出現數值問題……總之，參數真心不好調。
- 模擬收斂速度大概是200個球不到一分鐘吧。
- 還沒有用到 @小侯飛氘 提到的元胞列表法，用的話應該可以加快模擬收斂速度吧。
以上。
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~騙贊分割線~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
算了，估計也沒人點贊，還是直接給程序吧……主程序：鏈接：http://pan.baidu.com/s/1pL4g18v 密碼：mzhscircle.m函數：鏈接：http://pan.baidu.com/s/1gfBix2b 密碼：x9xw本來想試試看各位的演算法的，就先自己擼了個暴力的做參考，結果秒出（不到0.1秒），就懶得試別的演算法了23333
r = 40;
[xData, yData, rData] = deal([]);
while r &> 3
x = rand * 600; y = rand * 600;
n = 0;
while any((x - xData) .^ 2 + (y - yData) .^ 2 &< (r + rData) .^ 2)
x = rand * 600; y = rand * 600;
n = n + 1;
if n &> 200, r = r / 1.2; n = 0; end
end
xData(end + 1) = x; yData(end + 1) = y; rData(end + 1) = r;
end
t = (-15 : 15)" / 15 * pi;
patch(xData + cos(t) * rData, yData + sin(t) * rData, "w")
axis image off
考慮了一個稍微有點不同的問題：已知圓的個數及每個圓的半徑，找出一組圓心坐標使得這些圓互不重疊且盡量緊貼。
思路：當兩個圓有重疊時，給它們一個很大的排斥力；當兩個圓不重疊時，給它們一個相對小的吸引力；然後最小化總勢能。通過設計引力勢能等勢線的形狀，能大致控制最終結果的輪廓（像圓或者像方之類）。例子：
要最後像個圓形的話——
要最後像個方形的話——
Python 代碼：
import numpy as np
from scipy.optimize import minimize
from scipy.spatial.distance import pdist
import matplotlib.pyplot as plt
import os
def obj_grad_round(x, radii, pair_index_to_index_pair, repulsion_scale, attraction_scale):
n = len(radii)
centroids = np.reshape(x, (n, 2))
circles_pdist_sqrd = pdist(centroids, "sqeuclidean")
obj = 0
grad = np.zeros((n, 2))
for pidx, dist_sqrd in enumerate(circles_pdist_sqrd):
i, j = pair_index_to_index_pair[pidx]
touching_dist_sqrd = (radii[i] + radii[j]) ** 2
grad_inc = 4 * (dist_sqrd - touching_dist_sqrd) * (centroids[i, :] - centroids[j, :])
scale = repulsion_scale if dist_sqrd &< touching_dist_sqrd else attraction_scale
obj += scale * (dist_sqrd - touching_dist_sqrd) ** 2
grad[i, :] += scale * grad_inc
grad[j, :] -= scale * grad_inc
return obj, grad.flatten()
def obj_grad_square(x, radii, pair_index_to_index_pair, repulsion_scale, attraction_scale):
n = len(radii)
centroids = np.reshape(x, (n, 2))
circles_pdist_sqrd = pdist(centroids, "sqeuclidean")
obj = 0
grad = np.zeros((n, 2))
for pidx, dist_sqrd in enumerate(circles_pdist_sqrd):
i, j = pair_index_to_index_pair[pidx]
touching_dist_sqrd = (radii[i] + radii[j]) ** 2
diff = centroids[i, :] - centroids[j, :]
if dist_sqrd &< touching_dist_sqrd:
# repulsion
obj += repulsion_scale * (dist_sqrd - touching_dist_sqrd) ** 2
grad_inc = 4 * (dist_sqrd - touching_dist_sqrd) * diff
grad[i, :] += repulsion_scale * grad_inc
grad[j, :] -= repulsion_scale * grad_inc
else:
# attraction
diff_sqrd = diff ** 2
if max(diff_sqrd) &>= touching_dist_sqrd:
# do not overlap after projection, needs attraction
if diff_sqrd[0] &>= diff_sqrd[1]:
obj_inc = (diff_sqrd[0] - touching_dist_sqrd) ** 2
grad_inc = np.array([4 * (diff_sqrd[0] - touching_dist_sqrd) * diff[0], 0])
else:
obj_inc = (diff_sqrd[1] - touching_dist_sqrd) ** 2
grad_inc = np.array([0, 4 * (diff_sqrd[1] - touching_dist_sqrd) * diff[1]])
obj += attraction_scale * obj_inc
grad[i, :] += attraction_scale * grad_inc
grad[j, :] -= attraction_scale * grad_inc
return obj, grad.flatten()
def find_layout(radii, shape="square", repulsion_scale=1e6, attraction_scale=1, attraction_decay=0.3, overlapping_tol=0., seed=None, verbose=False):
assert shape == "round" or shape == "square"
obj_grad = obj_grad_round if shape == "round" else obj_grad_square
n = len(radii)
n_pairs = int(n * (n-1) / 2)
pair_index_to_index_pair = [None] * n_pairs
pidx = 0
for i in range(n-1):
for j in range(i+1, n):
pair_index_to_index_pair[pidx] = i, j
pidx += 1
if seed:
np.random.seed(seed)
if verbose:
options = {"disp": True}
history_folder = "history_%s" % shape
if not os.path.exists(history_folder):
os.mkdir(history_folder)
else:
options = None
overlapped = True
n_trials = 1
x0 = np.random.randn(2 * n)
while overlapped:
if verbose:
print("repulsion scale = %g attraction scale = %g" % (repulsion_scale, attraction_scale))
opt_result = minimize(obj_grad, x0, args=(radii, pair_index_to_index_pair, repulsion_scale, attraction_scale), jac=True, options=options)
centroids = np.reshape(opt_result.x, (n, 2))
if verbose:
if not opt_result.success:
print(opt_result.message)
draw_layout(centroids, radii, "%s/%02d_attraction=%g.png" % (history_folder, n_trials, attraction_scale))
np.savetxt("%s/%02d_attraction=%g.dat" % (history_folder, n_trials, attraction_scale), centroids, fmt="%g")
overlapped = check_overlapping(centroids, radii, overlapping_tol, verbose)
if overlapped:
attraction_scale *= attraction_decay
n_trials += 1
x0 = opt_result.x
else:
overlapped = False
return centroids
def check_overlapping(centroids, radii, tol, verbose=False):
n = len(radii)
for i in range(n-1):
for j in range(i+1, n):
dij = np.linalg.norm(centroids[i, :] - centroids[j, :])
rhs = radii[i] + radii[j] - tol
if dij &< rhs:
if verbose:
print("Circle {i} and circle {j} overlap ({dij} = dist(C_{i}, C_{j}) &< r_{i} + r_{j} - tol = {rhs}).".format(i=i, j=j, dij=dij, rhs=rhs))
return True
return False
def draw_layout(centroids, radii, not_show_and_save_to=None):
fig, ax = plt.subplots()
for centroid, radius in zip(centroids, radii):
circle = plt.Circle(centroid, radius)
ax.add_artist(circle)
x = centroids[:, 0]
y = centroids[:, 1]
xmin = min(x - radii)
xmax = max(x + radii)
ymin = min(y - radii)
ymax = max(y + radii)
plt.axis([xmin, xmax, ymin, ymax])
ax.set_aspect("equal")
if not_show_and_save_to:
plt.savefig(not_show_and_save_to)
else:
plt.show()
plt.close()
if __name__ == "__main__":
seed = 1
np.random.seed(seed)
n_circles = 100
radii = np.random.rand(n_circles) / 10 + 0.1 # [0.1, 0.2)
overlapping_tol = 1e-4
shape = "square" # "round" or "square"
verbose = False
layout = find_layout(radii, shape, overlapping_tol=overlapping_tol, verbose=verbose)
draw_layout(layout, radii, "final_%s.png" % shape)
np.savetxt("radii.dat", radii, fmt="%g")
np.savetxt("final_%s.dat" % shape, layout, fmt="%g")
.
grasshopper狗怒答
雖然有點取巧 不過感謝@孔德雨 的提示 新版的grasshopper確實有voronoi 的botton 不過grasshopper生成圓的原理暫時沒有根據voronoi vertices生成 但是注意到 有一個叫incircle的命令可以generate incircles in a triangle 需要的是三角形三個邊界點即可 這樣有兩種思路 一個是利用voronoi分解成幾個不重合的三角形 細分網格得到不重合的圓 這個具體我還沒有探索過....
第二個思路是利用delaunay網格直接生成三角形網格得到一系列不重合圓 經過一系列探索orz
得到如下結果
這個截圖是圓生成的base 即那個delaunay網格
不過我對第一種解法也很感興趣 有空探索一下
附上grasshopper電池圖
不知道為什麼會有runtime error 但是結果是沒錯的
以上
------------------------------一個小時候的分割線-----------------------------------------
感謝 @王喬志 的建議.....在每個vertex上塞了一個圓
原理是找到每個vertex 在它周圍的圓上project上的最短距離 以它為半徑生成圓保證不會交疊
這個說起來容易 做起來發現數據樹那是相當麻煩 以至於我一會兒graft一會flatten..... 以下是結果
這個是加上delaunay網格可以看到 每個三角形內生成的圓和每個vertex生成的圓互相不會交疊
然後友情附上grasshopper截圖....我連電池的習慣不太好 有時候會有點累贅 沒有coding基礎見諒
以上
------------------------更新的分割線--------------------------------
我來更新了！！！在過去的一年了學了kangaroo 發現有個新的解法 是利用kangaroo的 sphere collision來保證可以容納足夠多的圓而不重合 因為kangaroo是利用物理原理也就是說類似在一個容器里放進很多小球 所以他們一定不會重合....
截圖就不放了 感興趣的可以自己去了解一下kangaroo這個強大的插件
感謝 @王坤 在評論中的提醒，這應該是一個 Diffusion-limited aggregation 分形。粒子隨機聚集時就會表現出這種性質。--------------------------------更新了代碼，並添加了一些注釋---------------------------------在參數化設計課上用grasshopper做過一個類似的作業：如何用不同半徑的圓盡量鋪滿一個平面。分享一下我當時的解法：
- 先隨機生成兩個外切圓C1，C2；
- 生成一個隨機半徑r，分別在圓心O1, O2以半徑(R1+r), (R2+r)畫圓，得到第三個外切圓C3的圓心（隨機取其中一個交點）；
- 對C2, C3（或C1, C3）重複以上過程。
核心思路就是這樣，具體實現的話還要考慮不能和已有的圓相交以及邊界問題。
import Rhino.Geometry as rg
import rhinoscriptsyntax as rs
import random
outCircles = [] #記錄所有生成的圓
tanCircle = [] #記錄所有圓的相切關係和生成失敗次數
# 檢測圓是否在界內
def withInBound(circle):
if not(hasBound): return True
circle = circle.ToNurbsCurve()
curveList = rg.Curve.CreateBooleanIntersection(circle,bound)
if (len(curveList) &> 0 ):
curve = curveList[0]
return curve.EpsilonEquals(circle,10**-10)
else:
return False
# 檢測是否和已有的圓相交
def notIntersect(circle1):
center1 = circle1.Center
r1 = circle1.Radius
for circle2 in outCircles:
center2 = circle2.Center
r2 = circle2.Radius
dist = rs.Distance(center1, center2)
if dist&<(r1+r2):
return False
return True
#生成和C1,C2相切的圓
def generateCircle(index):
circle1, circle2, fails = tanCircle[index]
centerPoint = (circle1.Center + circle2.Center)/2
r = 1.5 * random.uniform(minRad, maxRad)
x = random.choice([0,1])
for i in range(3):
circle3 = rs.AddCircle(circle1.Center,circle1.Radius+r)
circle4 = rs.AddCircle(circle2.Center,circle2.Radius+r)
result = rs.CurveCurveIntersection(circle3,circle4) #求C3，C4交點（兩個）
rs.DeleteObjects([circle3, circle4])
if (result == None): continue
centerPoint = rg.Point3d(result[x][1])
circle5 = rg.Circle(centerPoint,r) #得到相切圓C5
if withInBound(circle5) and notIntersect(circle5):
tanCircle.append((circle1, circle5, 0)) #記錄相切關係
tanCircle.append((circle2, circle5, 0))
outCircles.append(circle5)
return True
r *= 0.8 #如不成功則嘗試縮小半徑
#r = 1.5 * random.uniform(minRad, maxRad)
x = not(x) #使用另外一個相切圓
tanCircle[index] = (circle1,circle2,fails+1)
return False
def debug():
print(i)
x = []
y = []
z = []
for circle in outCircles:
center = circle.Center
x += [center[0]]
y += [center[1]]
z += [center[2]]
fails= []
for circles in tanCircle:
fails += [(outCircles.index(circles[0]),outCircles.index(circles[1]),circles[2])]
return fails,x,y,z
if (hasBound): bound = rs.coercecurve(bound)
#第一個圓
center = rg.Point3d(0,0,0)
r = 1.5 * random.uniform(minRad, maxRad)
circle1 = rg.Circle(center,r)
outCircles.append(circle1)
#第二個圓
r = 1.5 * random.uniform(minRad, maxRad)
center = rg.Point3d(circle1.Radius+r,0,0)
circle2 = rg.Circle(center,r)
outCircles.append(circle2)
tanCircle.append((circle1,circle2,0))
#開始生成
for i in range(iteration):
for j in range(len(tanCircle)):
if (tanCircle[j][2] &< fails): generateCircle(j)
fails,x,y,z = debug()
問題本質：
首先這是一個不等原packing問題，即利用不等大的圓去填充給定的空間（給定的空間可以是四邊形，三角形或者圓形）。搜索關鍵字 「不等圓packing」可以找到相關資料。 「不等原packing」 是一個NP 難度問題。問題應用：
最簡單能想到的就是裁縫如何最大效率利用布料。解決方法：總結一下可能的解決方法，解決的問題和方式不完全和題目給定的一致，（提問者是填充矩形，但下面以填充圓形為例，填充給定大小和個數的圓）：1. 最大穴位演算法（簡而言之將圓與其他圓相切的方式以此放下去）
- 如何跳出填充中次優的解（有些小圓無法塞進去）
- 如何加速收斂因碰撞而產生的彈性勢能
[2] 基於格局變換策略的不等圓 Packing 問題求解演算法，計算機應用研究 [3] 解不等圓 packing 問題擬人擬物演算法初態選取，華中理工大學學報 [4] 支持求解原型 packing 問題的兩個擬人策略，中國科學 [5] 解 packing 及 CNF-SAT 問題的擬物擬人方法，華中理工大學學報 [6] 求解不等圓 packing 問題的一個啟發式演算法，計算機研究與發展
這問題前些年我其實考慮過，是在任意圖形中畫不相交的圓，思路大概如下：
- 得到一個任意形狀的區域
- 區域內任意生成一個點作圓心c
- 測量這個點到區域的最近距離d，並生成一個[0,d]值作為這個小圓的半徑r
- 從區域中摳掉這個小圓
- 重複2-4
當然，如果真的生成[0,d]的隨機半徑會發現遍地的極小圓，所以做了一丁點的修改，代碼如下
disk = Reap[
region =
RegionUnion[
BoundaryDiscretizeGraphics[
CountryData[#, "Polygon"]] /@ {"China", "Taiwan"}];
Do[p = RandomPoint[region];
rad = If[(tem = Abs[SignedRegionDistance[region, p]]) &< .2, tem,
RandomReal[{.2,
Min[{tem, Min@(Subtract @@ RegionBounds@region)/40}]}]];
region =
RegionDifference[region, DiscretizeRegion@Sow[Disk[p, rad]]],
2500]][[-1, -1]]; Graphics[
Transpose[{RandomColor[
Hue[1/3, NormalDistribution[.6, .2], NormalDistribution[.6, .07]],
disk // Length], disk}]]
我的暴力破解思路：
循環隨機生成一些圓的中心點，看看如果沒碰撞到其他圓的話，就在本次循環中長大一點點代碼：import javax.swing.*;
import java.awt.*;
import java.util.*;
public class ProduceVariableCircle extends JFrame{
public class Circle{
int x;//泡泡X坐標
int y;//泡泡y坐標
int r = 0;//泡泡半徑
double iSpeed;//泡泡膨脹速度
boolean bGrow = true;//泡泡還有膨脹空間嗎
}
MyPanel mp=null;
ArrayList&
public static void main(String[] args) {
ProduceVariableCircle pvc=new ProduceVariableCircle();
}
public ProduceVariableCircle(){
//循環2000次
for(int iRepeatTime=0;iRepeatTime&<2000;iRepeatTime++){
//每次隨機生成100個泡泡中心點，加到現存泡泡集合里
for(int iNewAdd=0;iNewAdd&<100;iNewAdd++){
Circle tempCircle = new Circle();
tempCircle.x = (int)(Math.random()*600);
tempCircle.y = (int)(Math.random()*600);
tempCircle.iSpeed = 1+(int)(Math.random()*20);
CircleList.add(tempCircle);
}
//逐一查看每個泡泡還有沒有膨脹空間
for(int i = CircleList.size()-1;i &>=0; i--){
//早就膨脹不了的，就不考慮了
if(!CircleList.get(i).bGrow){
continue;
}
//剩下的先膨脹一點點試試看
int r;
if(CircleList.get(i).r == 0){
r = (int)(CircleList.get(i).iSpeed);
}else{
r = (int)(CircleList.get(i).r*CircleList.get(i).iSpeed);
}
//一膨脹就馬上碰壁的，那就算了
int x = CircleList.get(i).x;
int y = CircleList.get(i).y;
if((x+r &> 600)||(x-r &< 0)||(y+r &> 600)||(y-r &< 0) ){
CircleList.get(i).bGrow = false;
continue;
}
//否則的話可以檢測和其他的泡泡有沒有撞到
for(int j = 0;j &< CircleList.size(); j ++){
//那些個根本就不在這個泡泡周圍的，應該不用考慮了吧（性能考慮）
int ox = CircleList.get(j).x;
int oy = CircleList.get(j).y;
int or = CircleList.get(j).r;
if((or&<=0)||(Math.abs(x-ox) &> r+ or)||(Math.abs(y-oy) &> r+ or)){
continue;
}
//在周圍的泡泡，那就看看撞到了沒有,撞到的話那就算了
double iDistence = Math.sqrt(Math.pow((x - ox),2)
+Math.pow((y - oy),2));
if((iDistence &< or)||(iDistence &< r + or)){
CircleList.get(i).bGrow = false;
break;
}
}
//一路過關斬將，終於發現確實是有膨脹空間的，那就快快長大吧
if(CircleList.get(i).bGrow){
CircleList.get(i).r = r;
}
//還是一個點，就沒有膨脹空間了的話，那隻好拜拜了
if(!CircleList.get(i).bGrow CircleList.get(i).r==0){
CircleList.remove(i);
}
}
}
mp =new MyPanel();
this.add(mp);
this.setSize(600, 600);
this.setLocation(600, 300);
this.setDefaultCloseOperation(this.EXIT_ON_CLOSE);
this.setVisible(true);
}
class MyPanel extends JPanel{
public void paint(Graphics g){
super.paint(g);
for(int i = 0;i &< CircleList.size(); i ++){
if(CircleList.get(i).r &>0 ){
g.drawOval((int)(CircleList.get(i).x-CircleList.get(i).r),
(int)(CircleList.get(i).y-CircleList.get(i).r),
(int)(2*CircleList.get(i).r),
(int)(2*CircleList.get(i).r));
}
}
}
}
}
上述代碼的效果如下：
我們以前考過類似的題……出題人應該是杜教。
可以這樣做：
- 我們先隨機生成一個遞減的序列，用來表示半徑。
- 接下來對於每個半徑，隨機出圓心的位置，如果能放下就直接放在那裡。這個隨機過程重複50次，如果每次都無法放下，則認為這個圓無法放下。
- 重複步驟1、2，直到無法操作為止。
gcj2012出過一道這樣的題，我當時做法是從左下角開始隨機圓心，然後做一個最大的外切圓，一直這樣直到平面滿了為止
且說說如何高效的判斷圓是否重疊。
很多物理模型中，高能粒子、原子、缺陷團簇乃至宏觀小球都採用圓或球來近似。當他們互相接近時會有相互作用，因此圓以及三維的球之間的碰撞判斷其實是非常有物理意義的一個問題。
前面很多答主給出了如何生成這些圓，但好像沒有答案討論如何高效的判斷圓與圓之間是否碰撞。而在分子動力學模擬中，有時需要對上百萬個原子進行上百萬步的模擬，每一步都要對每個原子進行這樣的判斷，判斷的效率自然是至關重要的。這裡我給出幾個演算法，拋磚引玉。
--------------------------------------------------
暴力枚舉法：每新添加一個圓，遍歷所有已存在的圓，判斷是否與新圓重疊。時間複雜度O(n)。
問題加了限定條件，此方法可能不再適用。====生成點，構造三角形，找內切圓。
唔....你們都欽點了邊界是矩形嗎....不必如此,雖然對於任意封閉曲線
positions=ImageValuePositions[EdgeDetect@img,White];
gra=Graphics[FilledCurve[Line[positions[[Last@FindShortestTour[positions]]]]]]
mesh=TriangulateMesh[BoundaryDiscretizeGraphics[gra,MaxCellMeasure-&>1],MaxCellMeasure-&>[Infinity]]
Graphics[MeshPrimitives[mesh,2]/.Polygon-&>Insphere]
imgBounds=Transpose[{{0,0},ImageDimensions[img]}];
vm=VoronoiMesh[ImageValuePositions[EdgeDetect[img,5],White],imgBounds];
vml=NestList[VoronoiMesh[Mean@@@MeshPrimitives[#,2],imgBounds],vm,3];
Graphics[Table[{RGBColor[ImageValue[img,Mean@@p]],p},{p,MeshPrimitives[Last[vml],2]}]]
那隻好找一些慢一點的演算法了,比如K近鄰.最大穴演算法什麼的....
唉...對呀...想歪了,其實沒有什麼邊界條件只有約束條件,求有界區域Distance=Compile[{{pt,_Real,1},{Centers,_Real,2},{Radii,_Real,1}},(Sqrt[Abs[First@#-First@pt]^2+Abs[Last@#-Last@pt]^2]/@Centers)-Radii,Parallelization-&>True,CompilationTarget-&>"C",RuntimeOptions-&>"Speed"];
shape=ColorNegate@EdgeDetect[img,2]
Centers=Radii={};
Module[{dim,dt,pt,m,d,r},dim=ImageDimensions[shape];
dt=DistanceTransform[shape];
TimeConstrained[While[True,While[While[pt=RandomReal[{1,#}]/@dim;
(m=3+ImageValue[dt,pt])&<5];
(d=Min[Distance[pt,Centers,Radii]]-2)&<5];
r=Min[10,d,m];Centers=Join[Centers,{pt}];
Radii=Join[Radii,{r}]],10]];
Graphics@MapThread[Disk,{Centers,Radii}]
Graphics[With[{col=RGBColor@ImageValue[img,RegionCentroid@#]},{EdgeForm@col,col,#}]/@MapThread[Disk,{Centers,Radii}]]
隨機生成點，算voronoi集，畫內切圓。
用一方碗 乘一碗帶油的湯，湯上面的圓都是不一樣大且不會重疊的！
圓推擠是個非常有效而且簡單的東西
哈，這個簡單！以下假定n為圓的數量。首先如何在規定的[上限,下限]半徑範圍內生成那麼多的不同半徑：只需用（上限-下限）*(k-1)/(n-1)+上限 即可，其中k為第幾個圓。若n為1那就xjb畫就好了。例如5個圓，半徑範圍為1-2mm，那麼則分別為1, 1.25, 1.5, 1.75, 2，保證不會有重複的半徑。做完半徑後，接下來可以生成圓了。首先，我們在一望無際的x軸上平鋪這些圓。
Mathematica的一種畫法，效率不高，主要是簡單
data=Table[{RandomReal[{0,10},2],RandomReal[{1/4,1}]},50];
disks=DeleteDuplicates[Disk@@@data,Not@*RegionDisjoint];
Graphics[{Circle@@@data,Opacity[0.2],disks}]
不邀自來~
%% Example
% This example illustrates how to draw circles randomly using MATLAB.
%% Preliminary work and settings
clc;
clear;
%Define the boundary of mian drawing domain
Xmax=1000;
Xmin=0;
Ymax=1000;
Ymin=0;
%The number of the circles needed to be drawn
Totalnum=200;
%The max and min radius values of a circle
Rmax=20;
Rmin=0;
%Initialization Parameters
Existnum=0; % The numbers of the existing circles
XPlist=[]; % The x position values of the center of existing circles
YPlist=[]; % The y position values of the center of existing circles
Rlist=[]; % The radius of the center of existing circles corresponding center positions
%% Drawing
figure;
box on;
axis equal;
hold on;
axis([Xmin Xmax Ymin Ymax]);
while 1
if Existnum&>Totalnum
break;
end
% Generate the radius and center position of a circle randomly
r=rand*(Rmax-Rmin)+Rmin;
xp=rand*(Xmax-Xmin)+Xmin;
yp=rand*(Ymax-Ymin)+Xmin;
% Check whether the standby circle is out of the boundarry
if xp+r&<=Xmax xp-r&>=Xmin yp+r&<=Ymax yp-r&>=Ymin
if Existnum==0
Existnum= Existnum+1;
rectangle("Position",[xp-r,yp-r,2*r,2*r],"Curvature",[1,1],"FaceColor",[rand,rand,rand]);
Rlist(Existnum)=r;
XPlist(Existnum)=xp;
YPlist(Existnum)=yp;
elseif Existnum~=0
% Check whether the standby circle is in existing circles
if all((XPlist-xp).^2+(YPlist-yp).^2&>=(Rlist+r).^2)
Existnum= Existnum+1;
rectangle("Position",[xp-r,yp-r,2*r,2*r],"Curvature",[1,1],"FaceColor",[rand,rand,rand]);
Rlist(Existnum)=r;
XPlist(Existnum)=xp;
YPlist(Existnum)=yp;
end
end
end
end
歡迎批評指正~
推薦閱讀：
※stl partition演算法有兩種寫法，哪種效率高？
※演算法導論第三版的翻譯水平如何？
※該直接上《演算法導論》 還是先看完 《演算法 第四版》？
※有能算出反三角函數的準確值的演算法嗎？
TAG:Python | 演算法 | MATLAB | WolframMathematica | Grasshopper |