Progressive Multijittered (0, 2) Sequence

pmg02.cpp

#include <iostream>
#include <random>

/* This implements the pseudo-code from http://graphics.pixar.com/library/ProgressiveMultiJitteredSampling/
 */

/* replace with your desired random number generator */
static std::random_device rd;
static std::mt19937 rng(rd());
static std::uniform_real_distribution<> dis(0, 1);
static float rnd() {
	float f;
	do {
		f = dis(rng);
	}
	while(!(f < 1.0f));
	return f;
}

typedef struct float2 {
	float x, y;
} float2;

class PMJ02_Generator {
public:
	static std::vector<float2> generate_2D(int size) {
		std::vector<float2> points;
		points.reserve(size);
		PMJ02_Generator g;
		float2 point;
		point.x = rnd();
		point.y = rnd();
		points.emplace_back(point);
		int N = 1;
		while (N < size) {
			g.extend_sequence_even(points, N);
			g.extend_sequence_odd(points, 2 * N);
			N = 4 * N;
		}
		return points;
	}

protected:
	PMJ02_Generator() {}

	void extend_sequence_even(std::vector<float2> &points, int N) {
		int n = sqrtf(N);
		mark_occupied_strata(points, N);
		for (int s = 0; s < N; ++s) {
			float2 oldpt = points[s];
			float i = floorf(n * oldpt.x);
			float j = floorf(n * oldpt.y);
			float xhalf = floorf(2.0f * (n * oldpt.x - i));
			float yhalf = floorf(2.0f * (n * oldpt.y - j));
			xhalf = 1.0f - xhalf;
			yhalf = 1.0f - yhalf;
			generate_sample_point(points, i, j, xhalf, yhalf, n, N);
		}
	}

	void extend_sequence_odd(std::vector<float2> &points, int N) {
		int n = sqrtf(N / 2);
		mark_occupied_strata(points, N);
		std::vector<float> xhalves(N / 2);
		std::vector<float> yhalves(N / 2);
		for (int s = 0; s < N / 2; ++s) {
			float2 oldpt = points[s];
			float i = floorf(n * oldpt.x);
			float j = floorf(n * oldpt.y);
			float xhalf = floorf(2.0f * (n * oldpt.x - i));
			float yhalf = floorf(2.0f * (n * oldpt.y - j));
			if (rnd() > 0.5f) {
				xhalf = 1.0f - xhalf;
			} else {
				yhalf = 1.0f - yhalf;
			}
			xhalves[s] = xhalf;
			yhalves[s] = yhalf;
			generate_sample_point(points, i, j, xhalf, yhalf, n, N);
		}
		for (int s = 0; s < N / 2; ++s) {
			float2 oldpt = points[s];
			float i = floorf(n * oldpt.x);
			float j = floorf(n * oldpt.y);
			float xhalf = 1.0f - xhalves[s];
			float yhalf = 1.0f - yhalves[s];
			generate_sample_point(points, i, j, xhalf, yhalf, n, N);
		}
	}

	void generate_sample_point(std::vector<float2> &points, float i, float j,
							   float xhalf, float yhalf, int n, int N) {
		int NN = 2 * N;
		float2 pt;
		do {
			pt.x = (i + 0.5f * (xhalf + rnd())) / n;
			pt.y = (j + 0.5f * (yhalf + rnd())) / n;
		} while (is_occupied(pt, NN));
		mark_occupied_strata1(pt, NN);
		points.emplace_back(pt);
	}

	void mark_occupied_strata(std::vector<float2> &points, int N) {
		int NN = 2 * N;
		int num_shapes = log2f(NN) + 1;
		occupiedStrata.resize(num_shapes);
		for (int shape = 0; shape < num_shapes; ++shape) {
			occupiedStrata[shape].resize(NN);
			for (int n = 0; n < NN; ++n) {
				occupiedStrata[shape][n] = false;
			}
		}
		for (int s = 0; s < N; ++s) {
			mark_occupied_strata1(points[s], NN);
		}
	}

	void mark_occupied_strata1(float2 pt, int NN) {
		int shape = 0;
		int xdivs = NN;
		int ydivs = 1;
		do {
			int xstratum = xdivs * pt.x;
			int ystratum = ydivs * pt.y;
			occupiedStrata[shape][ystratum * xdivs + xstratum] = true;
			shape = shape + 1;
			xdivs = xdivs / 2;
			ydivs = ydivs * 2;
		} while (xdivs > 0);
	}

	bool is_occupied(float2 pt, int NN) {
		int shape = 0;
		int xdivs = NN;
		int ydivs = 1;
		do {
			int xstratum = xdivs * pt.x;
			int ystratum = ydivs * pt.y;
			if (occupiedStrata[shape][ystratum * xdivs + xstratum]) {
				return true;
			}
			shape = shape + 1;
			xdivs = xdivs / 2;
			ydivs = ydivs * 2;
		} while (xdivs > 0);
		return false;
	}

private:
	std::vector<std::vector<bool>> occupiedStrata;
};

int main(int argc, char **argv) {
	constexpr int num_points = 1024;
	std::vector<float2> points = PMJ02_Generator::generate_2D(num_points);
	for (auto point : points) {
		std::cout << point.x << ", " << point.y << "\n";
	}
}