Generowanie próbek tą metodą jest wyjątkowo wolne. Złożoność obliczeniowa to $N log N$ W kodzie używam i tak już zmodyfikowanej metody, której złożoność obliczeniowa to N.
Rozkład punktów wygenerowanych przez sampler jest porównywalny do rozkładu fotoreceptorów siatkówki.
Aby wygenerować zbiór punktów o rozkładzie Poisson-Disc bierzemy losowy punk i kładziemy go w przestrzeni którą chcemy zapełnić (może być wielowymiarowa). Następnie losujemy pewien zbiór punktów (liczba prób zakończona w kodzie jako
Repeats) i wybieramy ten, który jest najbardziej oddalony od już postawionych (minimalny dystans do już postawionych jest maksymalny) oraz nie jest bliżej niż pewien promień $r_{MIN}$, który w praktyce wyznacza nam gęstość rozkładu. Algorytm ten jest bardzo wolny. Po skończonym czasie może się zdarzyć (powinno się zdarzyć jeśli mamy na myśli skończony bardzo długi czasu, że nie możemy już dodać żadnych nowych punktów. W praktyce algorytm trzeba przerwać dużo wcześniej, co z drugiej strony może dać nam rozkład z dziurami. Możemy sobie wyobrazić metodę w której losowanie odbywa się tylko w obszarach, w których mamy pewność, że dziury są. W szybszej wersji tworzymy siatkę o rozmiarze tak dobranym by średnio jeden punkt wpadał w jedno oczko siatki i w przypadku wylosowania punktu sprawdzamy tylko obszar 3x3 w takiej siatce. Rozmiar oczka siatki wyznaczany jest przez $r_{MIN}$. Dwa punkty w sąsiednich oczkach siatki (także po przekątnych) maksymalnie od siebie oddalone nie zostawią w środku miejsca na kolejny punkt. Ponieważ odległość taka jest największa pod przekątnej rozmiar siatki definiuje się jako: $\frac{r_{MIN}}{\sqrt{2}}$.
Na początku losujemy punkt i dodajemy go do siatki i do listy aktywnej. Następnie tak długo jak lista aktywna nie jest pusta losujemy z niej jeden punkt A. Dla tego punktu losujemy
Repeats punktów C wokół niego oddalonych od $0.5 \cdot r_{MIN}$ do $r_{MIN}$. Następnie dla każdego punktu B z C sprawdzamy siatkę 3x3 wokół punktu B. Jeśli punkt B jest oddalony od wszystkich punktów w siatce 3x3 o więcej niż $r_{MIN}$ wstawiamy go siatki i do listy aktywnej. W przypadku gdy nie udało nam się dodać żadnego nowego punktu wokół A usuwamy go z listy aktywnej. Tym samym czym większa liczba prób tym mniej potencjalnych dziur, które da się jeszcze coś wepchnąć. Przykłady rozkładu próbek:
Kod samplera Poisson:
public class PoissonSampler : NonUniformSampler
{
public int Repeats = 20;
public double Density = 1;
private double m_cell_size;
private double m_radius;
private int m_grid_width;
private int m_grid_height;
private Random m_random;
private int m_width;
private int m_height;
private int[,] m_grid;
List<int> m_actives = new List<int>();
List<Vector2> m_samples = new List<Vector2>();
private void PrepareSamples()
{
m_width = Film.Width;
m_height = Film.Height;
if (Film.Scene.RenderOptions.Renderer is SamplesVisualizerRenderer)
{
double scale = (Film.Scene.RenderOptions.Renderer as
SamplesVisualizerRenderer).Scale;
m_width = (m_width / scale).Round();
m_height = (m_height / scale).Round();
}
m_radius = 1.0 / Density;
m_cell_size = m_radius / Math.Sqrt(2);
m_grid_width = (int)(m_width / m_cell_size) + 1;
m_grid_height = (int)(m_height / m_cell_size) + 1;
m_grid = new int[m_grid_width, m_grid_height];
m_grid.Fill(-1);
AddFirstPoint();
while (m_actives.Count != 0)
{
int index = m_random.Next(m_actives.Count);
bool found = false;
for (int i = 0; i < Repeats; i++)
found |= AddNextPoint(m_actives[index]);
if (!found)
m_actives.RemoveAt(index);
}
}
private Vector2 GenerateRandomAround(int a_index)
{
double radius = (m_radius + m_radius * m_random.NextDouble());
double angle = 2 * MathExtensions.PI * m_random.NextDouble();
Vector2 center = m_samples[a_index];
return new Vector2(center.X + radius * Math.Sin(angle),
center.Y + radius * Math.Cos(angle));
}
private bool AddNextPoint(int a_index)
{
Vector2 q = GenerateRandomAround(a_index);
if ((q.X < 0) || (q.X >= m_width) || (q.Y < 0) || (q.Y >= m_height))
return false;
int index_x = GetXIndex(q);
int index_y = GetYIndex(q);
for (int x = Math.Max(0, index_x - 2);
x < Math.Min(m_grid_width, index_x + 3); x++)
{
for (int y = Math.Max(0, index_y - 2);
y < Math.Min(m_grid_height, index_y + 3); y++)
{
int index = m_grid[x, y];
if (index != -1)
{
if ((q - m_samples[index]).Length < m_radius)
return false;
}
}
}
m_grid[index_x, index_y] = m_samples.Count;
m_actives.Add(m_samples.Count);
m_samples.Add(q);
return true;
}
private int GetXIndex(Vector2 a_point)
{
return (int)(a_point.X / m_cell_size);
}
private int GetYIndex(Vector2 a_point)
{
return (int)(a_point.Y / m_cell_size);
}
private void AddFirstPoint()
{
Vector2 p = new Vector2(m_width * m_random.NextDouble(),
m_height * m_random.NextDouble());
m_grid[GetXIndex(p), GetYIndex(p)] = m_samples.Count;
m_actives.Add(m_samples.Count);
m_samples.Add(p);
}
public override IEnumerable<Vector2> GetSamples(Rectangle a_rect)
{
int top = (int)(a_rect.Top / m_cell_size) + 1;
int bottom = (int)(a_rect.Bottom / m_cell_size) + 1;
int left = (int)(a_rect.Left / m_cell_size) + 1;
int right = (int)(a_rect.Right / m_cell_size) + 1;
for (int x = Math.Max(0, left - 2);
x <= Math.Min(m_grid_width - 1, right + 2); x++)
{
for (int y = Math.Max(0, top - 2);
y <= Math.Min(m_grid_height - 1, bottom + 2); y++)
{
if (m_grid[x, y] != -1)
{
Vector2 p = m_samples[m_grid[x, y]];
if (p.X < a_rect.Left)
continue;
if (p.X > a_rect.Right)
continue;
if (p.Y < a_rect.Top)
continue;
if (p.Y > a_rect.Bottom)
continue;
yield return p;
}
}
}
}
}
Parametr skalowania można pominąć, jest on używany tylko podczas wizualizacji rozkładu próbek i w periodogramie.
Brak komentarzy:
Prześlij komentarz