/*
								+----------------------------------+
								|                                  |
								| ***  Poisson disc generator  *** |
								|                                  |
								|   Copyright © -tHE SWINe- 2007   |
								|                                  |
								|           Poisson.cpp            |
								|                                  |
								+----------------------------------+
*/

/**
 *	@file lml/Poisson.cpp
 *	@brief Poisson disc generator
 *	@date 2007
 *	@author -tHE SWINe-
 */

#include "../NewFix.h"
#include "../CallStack.h"
#include <vector>
#include <algorithm>
#include <functional>
#include <math.h>
#include "../MinMax.h"
#include "Poisson.h"

#if defined(_MSC_VER) && !defined(__MWERKS__) && !defined(for) && _MSC_VER <= 1200
#define for if(0) {} else for
#endif // _MSC_VER && !__MWERKS__ && !for && _MSC_VER <= 1200
// msvc 'for' scoping hack

/*
 *								=== TDisc ===
 */

TDisc::TDisc(const Vector2f &r_v_center, float _f_radius)
	:v_center(r_v_center), f_radius(_f_radius), n_interval_num(1)
{
	p_arc_angle[0] = TInterval<float>(0, float(2 * f_pi));
}

bool TDisc::Subtract(const TDisc &r_t_disc)
{
	_ASSERTE(r_t_disc.f_radius == f_radius);

	Vector2f v_diff = r_t_disc.v_center - v_center;
	float f_len = v_diff.f_Length();

	if(f_len >= 2 * f_radius)
		return false;
	// discs miss each other

	_ASSERTE(f_len > 0);
	// do not handle situation when discs exactly overlap

	float f_alpha = (float)acos(f_len / (2 * f_radius));
	float f_angle = float(atan2(v_diff.y, v_diff.x));
	// calculate incident angle and intersection disc (angular) length

	if(f_angle < f_alpha)
		f_angle += float(2 * f_pi);

	float f_angle0 = f_angle - f_alpha;
	float f_angle1 = f_angle + f_alpha;
	// calculate angle interval

	SubtractAngleInterval(f_angle0, f_angle1);
	// subtract angular range from this disc

	return true;
}

void TDisc::SubtractAngleInterval(float f_angle0, float f_angle1)
{
	_ASSERTE(f_angle0 >= 0 && f_angle0 <= 2 * f_pi);
	_ASSERTE(f_angle1 >= f_angle0);
	_ASSERTE(f_angle1 >= 0 && f_angle1 <= 4 * f_pi);
	// make sure angles are in "nice" intervals

	for(;;) {
		float f_max = min(f_angle1, float(2 * f_pi));
		TInterval<float> t_sub_interval(f_angle0, f_max);

		TInterval<float> p_new_arcs[n_max_interval_num];
		int n_arc_num = 0;

		for(int i = 0, n = n_interval_num; i < n; ++ i) {
			if(p_arc_angle[i].b_Disjunct(t_sub_interval)) {
				p_new_arcs[n_arc_num] = p_arc_angle[i];
				if(++ n_arc_num == n_max_interval_num)
					break;
				continue;
			}
			// intervals do not intersect

			if(p_arc_angle[i].f_min < t_sub_interval.f_min) { // not equal, then it could be empty
				_ASSERTE(!p_arc_angle[i].t_LeftComplement(t_sub_interval).b_Empty());
				p_new_arcs[n_arc_num] = p_arc_angle[i].t_LeftComplement(t_sub_interval);
				if(++ n_arc_num == n_max_interval_num)
					break;
			}
			if(p_arc_angle[i].f_max > t_sub_interval.f_max) { // not equal, then it could be empty
				_ASSERTE(!p_arc_angle[i].t_RightComplement(t_sub_interval).b_Empty());
				p_new_arcs[n_arc_num] = p_arc_angle[i].t_RightComplement(t_sub_interval);
				if(++ n_arc_num == n_max_interval_num)
					break;
			}
			// subtract intervals (possible subdivision of interval to up to two intervals)
		}
		// subtract intervals

		n_interval_num = n_arc_num;
		memcpy(p_arc_angle, p_new_arcs, sizeof(TInterval<float>) * n_arc_num);
		// copy intervals back (presumably cheaper than redistributing)

		if(f_angle1 <= 2 * f_pi)
			break;
		f_angle0 = 0;
		f_angle1 -= float(2 * f_pi);
	}
	// subtract intervals, subdivide interval going over 2pi
}

float TDisc::f_FreeArc_Length() const
{
	float f_result = 0;
	for(int i = 0, n = n_interval_num; i < n; ++ i) {
		_ASSERTE(!p_arc_angle[i].b_Empty());
		f_result += p_arc_angle[i].f_max - p_arc_angle[i].f_min;
	}
	return f_result;
}

Vector2f TDisc::v_FreeArc_Point(float f_angle) const
{
	for(int i = 0, n = n_interval_num; i < n; ++ i) {
		_ASSERTE(!p_arc_angle[i].b_Empty());
		float f_arc_len = p_arc_angle[i].f_max - p_arc_angle[i].f_min;
		if(f_angle < f_arc_len) {
			_ASSERTE(f_angle >= 0);
			f_angle += p_arc_angle[i].f_min;
			return v_center + Vector2f(cos(f_angle), sin(f_angle)) * f_radius;
		}
		f_angle -= f_arc_len;
	}
	return v_center; // return center to mark error
}

/*
 *								=== ~TDisc ===
 */

/*
 *								=== CPoissonDisc_Set ===
 */

class CSumArcLengths {
protected:
	float m_f_sum;

public:
	inline CSumArcLengths()
		:m_f_sum(0)
	{}

	inline void operator ()(const TDisc &r_t_disc)
	{
		m_f_sum += r_t_disc.f_FreeArc_Length();
	}

	inline operator float() const
	{
		return m_f_sum;
	}
};

class CFindPointOnArc {
protected:
	float &m_r_f_arc_angle;

public:
	inline CFindPointOnArc(float &r_f_arc_angle)
		:m_r_f_arc_angle(r_f_arc_angle)
	{}

	inline bool operator ()(const TDisc &r_t_disc)
	{
		float f_arc_len = r_t_disc.f_FreeArc_Length();
		if(m_r_f_arc_angle < f_arc_len)
			return true;
		m_r_f_arc_angle -= f_arc_len;
		return false;
	}
};

CPoissonDisc_Set::CPoissonDisc_Set(float f_disc_radius,
	bool b_repeatable, const CRandomGeneratorModel &r_generator)
	:m_b_status(true)
{
	try {
		std::vector<TDisc> mirror_disc_list;
		// side-mirrored discs to assure repeatability of generated set

		std::vector<TDisc> spent_disc_list;
		// list of discs that do not own any arc length

		for(bool b_first_pass = true;; b_first_pass = false) {
			Vector2f v_new_center;
			if(b_first_pass) {
				v_new_center = Vector2f(r_generator.f_Rand0(), r_generator.f_Rand0());
				// randomly generate position of the first disc
			} else {
				float f_free_arc_sum = std::for_each(m_disc_list.begin(),
					m_disc_list.end(), CSumArcLengths());
				if(f_free_arc_sum < 1e-3f) {
					m_disc_list.insert(m_disc_list.end(),
						spent_disc_list.begin(), spent_disc_list.end());
					// put the other discs back

					//m_disc_list = mirror_disc_list; // debug

					return;
				}
				// there's no free arc to be occupied

				float f_free_arc_pos = f_free_arc_sum * r_generator.f_Rand0();
				// get angular position on free arc where new point should be generated

				std::vector<TDisc>::iterator p_disc_iter = std::find_if(m_disc_list.begin(),
					m_disc_list.end(), CFindPointOnArc(f_free_arc_pos));
				if(p_disc_iter == m_disc_list.end()) { // note that this could be due to a numerical error
					m_b_status = false;
					return;
				}
				// find the disc, f_free_arc_pos is modified so it now points to position
				// on the disc

				_ASSERTE(f_free_arc_pos >= 0 && f_free_arc_pos <=
					(*p_disc_iter).f_FreeArc_Length());
				v_new_center = (*p_disc_iter).v_FreeArc_Point(f_free_arc_pos);
				// find a new point on the disc
			}

			TDisc t_new_disc(v_new_center, f_disc_radius);
			// new disc

			const size_t n_mirror_disc_num_before = mirror_disc_list.size();
			if(b_repeatable) {
				for(int n_dy = -1; n_dy < 2; ++ n_dy) {
					float f_y = v_new_center.y + n_dy;
					if((f_y <= -f_disc_radius || f_y >= 1 + f_disc_radius))
						continue;
					for(int n_dx = -1; n_dx < 2; ++ n_dx) {
						float f_x = v_new_center.x + n_dx;
						if((n_dx != 0 || n_dy != 0) &&
						   (f_x > -f_disc_radius && f_x < 1 + f_disc_radius))
							mirror_disc_list.push_back(TDisc(Vector2f(f_x, f_y), f_disc_radius));
					}
				}
				// create several more discs to subtract from existing
				// discs to mimic wrap-arround behavior
			}

			if(v_new_center.x < f_disc_radius) {
				t_new_disc.Subtract(TDisc(Vector2f(-v_new_center.x,
					v_new_center.y), f_disc_radius));
			}
			if(v_new_center.x > 1 - f_disc_radius) {
				t_new_disc.Subtract(TDisc(Vector2f(2 - v_new_center.x,
					v_new_center.y), f_disc_radius));
			}
			if(v_new_center.y < f_disc_radius) {
				t_new_disc.Subtract(TDisc(Vector2f(v_new_center.x,
					-v_new_center.y), f_disc_radius));
			}
			if(v_new_center.y > 1 - f_disc_radius) {
				t_new_disc.Subtract(TDisc(Vector2f(v_new_center.x,
					2 - v_new_center.y), f_disc_radius));
			}
			// subtract borders (do it a better way?)
			// note that the branches are pairwise-exclusive for discs with radius less than 0.5

			for(size_t j = 0, n = mirror_disc_list.size(); j < n; ++ j)
				t_new_disc.Subtract(mirror_disc_list[j]);
			for(size_t j = 0, n = spent_disc_list.size(); j < n; ++ j)
				t_new_disc.Subtract(spent_disc_list[j]);
			// subtract mirrors and spent discs from it

			for(size_t i = 0, n = m_disc_list.size(); i < n; ++ i) {
				t_new_disc.Subtract(m_disc_list[i]);
				// from the new disc

				for(size_t j = n_mirror_disc_num_before, m = mirror_disc_list.size(); j < m; ++ j)
					m_disc_list[i].Subtract(mirror_disc_list[j]);
				bool b_subtracted = m_disc_list[i].Subtract(t_new_disc);
				// from other active discs

				if(b_subtracted && m_disc_list[i].f_FreeArc_Length() < 1e-3f) {
					spent_disc_list.push_back(m_disc_list[i]);
					m_disc_list.erase(m_disc_list.begin() + i);
					-- i;
					-- n;
				}
				// put spent discs away (saves subtracting from them)
			}
			// subtract discs from each other (O(n2), could accelerate by uniform grid cache)

			if(t_new_disc.f_FreeArc_Length() > 1e-3f)
				m_disc_list.push_back(t_new_disc);
			else
				spent_disc_list.push_back(t_new_disc);
			// add new disc
		}
	} catch(std::bad_alloc&) {
		m_b_status = false;
		// not enough memory
	}
}


/*
 *								=== ~CPoissonDisc_Set ===
 */

/*
 *								=== CSamplePattern ===
 */

CSamplePattern::CSamplePattern(size_t n_point_num, size_t n_point_num_tolerance,
	bool b_repeatable, CRandomGeneratorModel &r_generator, int n_seed)
{
	GenerateDiscs(n_point_num, n_point_num_tolerance, b_repeatable, r_generator, n_seed);
}

void CSamplePattern::Center()
{
	if(m_sample_list.empty())
		return;

	Vector2f v_min(m_sample_list[0]), v_max(m_sample_list[0]);
	for(size_t i = 0, n = m_sample_list.size(); i < n; ++ i) {
		for(int j = 0; j < 2; ++ j) {
			float f = m_sample_list[i][j];
			v_min[j] = (v_min[j] < f)? v_min[j] : f;
			v_max[j] = (v_max[j] > f)? v_max[j] : f;
		}
	}
	// determine bounding box

	Vector2f v_shift = (v_min + v_max) * -.5f;
	// determine shift

	for(size_t i = 0, n = m_sample_list.size(); i < n; ++ i)
		m_sample_list[i] += v_shift;
	// shift samples
}

void CSamplePattern::Scale(float f_scale)
{
	for(int i = 0, n = m_sample_list.size(); i < n; ++ i)
		m_sample_list[i] *= f_scale;
	// scale samples
}

bool CSamplePattern::GenerateDiscs(size_t n_point_num, size_t n_point_num_tolerance,
	bool b_repeatable, CRandomGeneratorModel &r_generator, int n_seed)
{
	float f_size = 1 / max(float(sqrt(double(n_point_num))), 1.0f);
	for(float f_step = f_size * .5f; f_step > 1e-10f; f_step *= .5f) {
		r_generator.Seed(n_seed);
		CPoissonDisc_Set pd_set(f_size, b_repeatable, r_generator);
		if(!pd_set.b_Status())
			return false;

		size_t n_cur_point_num = pd_set.r_DiscList().size();
		// generate set of points

		if(max(n_cur_point_num, n_point_num) - min(n_cur_point_num,
		   n_point_num) <= n_point_num_tolerance) {
			const std::vector<TDisc> &r_disc_list = pd_set.r_DiscList();
			std::insert_iterator<std::vector<Vector2f> >
				p_insert_iter(m_sample_list, m_sample_list.begin());
			std::copy(r_disc_list.begin(), r_disc_list.end(), p_insert_iter);
			// copy & convert

			return true;
		}
		else if(n_cur_point_num > n_point_num)
			f_size += f_step;
		else
			f_size -= f_step;
		// decide next action
	}
	// try to generate set with required number of points by altering disc radius

	/*r_generator.Seed(n_seed);
	for(int i = 0; i < 100; ++ i) {
		CPoissonDisc_Set pd_set(f_size, r_generator);
		if(!pd_set.Get(r_set))
			return false;
		int n_cur_point_num = r_set.size();
		// generate set of points using the same generator
		// several times (effectively changing seed), try luck

		if(abs(n_cur_point_num - n_point_num) <= n_tolerance)
			return true;
	}*/
	// try luck in case it's hopeless (can't balance
	// number of points just by setting disc radius)

	return false;
}

/*
 *								=== ~CSamplePattern ===
 */
