/* -*- C++ -*-
 * Copyright ©2004 Hugo Mills <hugo@carfax.org.uk>
 *
 * This software is distributed under the terms of the GNU GPL v3
 * For more information on the GPL, see the file COPYING or 
 * visit http://www.gnu.org/
 *
 * This software is distributed without warranty
 */

#include "display_geodesic.h"

#include <iostream>

void Display_Geodesic::plot(void)
{
	dvec p0(3);
	dvec p1(3);

	if(sphere->f(p0, pt0)
	   && sphere->f(p1, pt1))
	{
		dvec shifted_pos(3);
		dvec map_pos(2);
		ivec vp_pos0(2);
		ivec vp_pos1(2);

		bool p0map = 
			orbit->f(shifted_pos, p0)
			&& projection->f(map_pos, shifted_pos)
			&& viewport->f(vp_pos0, map_pos);
		bool p1map = 
			orbit->f(shifted_pos, p1)
			&& projection->f(map_pos, shifted_pos)
			&& viewport->f(vp_pos1, map_pos);

		midpoint_div(p0, p1, 
					 p0map, vp_pos0,
					 p1map, vp_pos1, 
					 0);
	}
}

void Display_Geodesic::midpoint_div(
	const dvec& p0, const dvec& p1,
	bool p0map, const ivec& vp_pos0, 
	bool p1map, const ivec& vp_pos1, int depth) const
{
	dvec shifted_pos(3);
	dvec map_pos(2);
	ivec vp_mid(2);
	bool midmap = true;
	bool h0 = true;
	bool h1 = true;

	dvec mid(3);
	split(mid, p0, p1);

	if(depth > 20)
		return;

	if(orbit->f(shifted_pos, mid)
	   && projection->f(map_pos, shifted_pos)
	   && viewport->f(vp_mid, map_pos))
	{
		if(p0map)
		{
			// 5 == arbitrary flatness, in pixels
			if(abs(vp_pos0[0] - vp_mid[0]) + abs(vp_pos0[1] - vp_mid[1]) <= 5)
			{
				bitmap->setpixel(vp_mid, colour);
				bitmap->setpixel(vp_pos0, colour);
				h0 = false;
			}
		}

		if(p1map)
		{
			if(abs(vp_pos1[0] - vp_mid[0]) + abs(vp_pos1[1] - vp_mid[1]) <= 5)
			{
				bitmap->setpixel(vp_mid, colour);
				bitmap->setpixel(vp_pos1, colour);
				h1 = false;
			}
		}
	}
	else
		midmap = false;

	if(h0) midpoint_div(p0, mid,
						p0map, vp_pos0,
						midmap, vp_mid, 
						depth+1);
	if(h1) midpoint_div(mid, p1, 
						midmap, vp_mid, 
						p1map, vp_pos1, 
						depth+1);
}

void Display_Geodesic::split(
	dvec& result, 
	const dvec& p0, 
	const dvec& p1) const
{
	result = p0 + p1;
	double r = sqrt(result[0]*result[0]
					+ result[1]*result[1]
					+ result[2]*result[2]);
	result /= r;
}