3

次の制約を使用して、正方形に N 個の点 (通常は 1e2 - 1e4) の 2D セットを作成したいと考えています。

  • すべてのポイント間に最小限の距離が必要です (ハードコア除外ゾーン)

  • 固定密度を取得したいので、正方形を埋めるポイントの数は事前に (または概算で) 与えられます (必要に応じて、後で正方形のサイズを少し調整できます)。

  • パターンは合理的に「ランダム」でなければなりません

  • 迅速な解決策が優先されます

以前はパッケージ spatstat で rStrauss を使用していましたが、特定のポイント数を確実に取得する方法を理解できず、おそらくタスクが難しすぎたため、関数がマシンを 10 分間停止させることがよくありました。これにはもっと適切な機能があるかもしれないと思います。

## regular grid of 1e2 points in [-10, 10]^2
xy = expand.grid(x=seq(-10, 10, length=10), y=seq(-10, 10, length=10))
N = NROW(xy)

編集：答えで示唆されているように

xyr = rSSI(r=0.1, N, win = owin(c(-10,10),c(-10,10)), N)
plot(xyr)

pp

4

2 に答える 2

3

rSSIもspatstatパッケージに含まれており、標準によっては速度を除いて、問題を処理します。それは筋金入りの抑制距離を持ち、目標数のポイントに到達します (または試行をあきらめます。あきらめるためのしきい値を設定できます)。配置はランダムです。特に速いとは思いませんが、約30秒で1e6阻害距離の単位正方形に点を作ることができました。1e-4速度と成功率は、ポイント数に対する阻害距離に大きく依存します。

于 2011-11-19T22:35:52.680 に答える
2

ほとんどの場合、Rcpp についてさらに学習する口実として、これを行うための小さな関数を試してみます。

require(inline)
require(Rcpp)

randPoints = cxxfunction(signature(r_n='int', r_mindist='float', r_maxiter='int'), body = 
' 
  using namespace std;

  int n = as<int> (r_n);
  float mindist = as<float> (r_mindist);
  int maxiter = as<int> (r_maxiter);

  RNGScope scope;
  bool tooclose;
  int iter;
  NumericVector rands (2);
  NumericMatrix points (n, 2);
  NumericVector dist (2);

  for (int i=0; i < n; i++) {
    iter = 0;
    do {
      iter++;
      tooclose = false;
      rands = runif(2, 0, 1);
      for (int idist=0; idist < i; idist++) {
        dist = rands - points(idist, _);
        dist = dist * dist;
        if (sqrt(accumulate(dist.begin(), dist.end(), 0.0)) < mindist) {
          tooclose = true;
          break;
        }
      }
    } while (tooclose & iter < maxiter);
    if (iter == maxiter) {
      Rprintf("%d iterations reached\\nOnly %d points found\\n", maxiter, i+1);
      break;
    }
    NumericMatrix::Row target(points, i);
    target = rands;
  }

  return(wrap(points));
'
, plugin='Rcpp')

次に、次のように使用できます。

> x = randPoints(1000, 0.05, 10000)
10000 iterations reached
Only 288 points found

そして、ここにプロットがあります：

x = x[as.logical(rowMeans(x != 0)), ]
dev.new(width=4, height=4)
plot(x)

ここに画像の説明を入力

于 2011-11-21T00:19:00.380 に答える