псевдослучайный узор с фиксированной плотностью и зоной отчуждения

Я хочу создать двумерный набор из N точек (обычно 1e2 - 1e4) в квадрате со следующими ограничениями:

  • между всеми точками должно быть минимальное расстояние (зона отчуждения жесткого сердечника)

  • количество точек, заполняющих квадрат, указано заранее (или приблизительная оценка), так как я хочу получить фиксированную плотность (при необходимости я могу потом немного изменить размер квадрата).

  • узор должен быть достаточно "случайным"

  • предпочтительнее быстрое решение

Раньше я использовал rStrauss в пакете spatstat, но я никогда не мог понять, как надежно получить заданное количество баллов, и довольно часто функция останавливала мою машину на 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


r spatial spatstat
person baptiste    schedule 19.11.2011    source источник

Ответы (2)


arrow_upward
3
arrow_downward

rSSI, также входящий в пакет spatstat, позаботится о ваших проблемах, кроме, возможно, скорости, в зависимости от ваших стандартов. Он имеет жесткую дистанцию ​​запрета и наберет целевое количество очков (или откажитесь от попыток - вы можете установить порог отказа), а места размещения являются случайными. Я не думаю, что это особенно быстро, но я смог создать 1e6 точек в единичном квадрате с расстоянием запрета 1e-4 примерно за 30 секунд. Скорость и вероятность успеха будут сильно зависеть от вашего расстояния торможения по отношению к количеству очков.

person Gregor Thomas    schedule 19.11.2011
comment
Интересная проблема. Мне было бы интересно узнать, какой алгоритм в целом является наиболее эффективным. Спасибо за ссылки! - person John Colby; 20.11.2011

arrow_upward
2
arrow_downward

В основном в качестве предлога узнать больше о 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)

введите описание изображения здесь

person John Colby    schedule 21.11.2011