#include <2geom/d2.h>
#include <2geom/sbasis.h>
#include <2geom/sbasis-geometric.h>
#include <2geom/orphan-code/intersection-by-smashing.h>

#include <cstdlib>
#include <cstdio>
#include <vector>
#include <algorithm>

namespace Geom {
using namespace Geom;

/*
 * Computes the top and bottom boundaries of the L_\infty neighborhood
 * of a curve. The curve is supposed to be a graph over the x-axis.
 */
static
void computeLinfinityNeighborhood( D2<SBasis > const &f, double tol, D2<Piecewise<SBasis> > &topside, D2<Piecewise<SBasis> > &botside ){
	double signx = ( f[X].at0() > f[X].at1() )? -1 : 1;
	double signy = ( f[Y].at0() > f[Y].at1() )? -1 : 1;

	Piecewise<D2<SBasis> > top, bot;
	top = Piecewise<D2<SBasis> > (f);
	top.cuts.insert( top.cuts.end(), 2);
	top.segs.insert( top.segs.end(), D2<SBasis>(SBasis(Linear( f[X].at1(), f[X].at1()+2*tol*signx)),
			                                    SBasis(Linear( f[Y].at1() )) ));
	bot = Piecewise<D2<SBasis> >(f);
	bot.cuts.insert( bot.cuts.begin(), - 1 );
	bot.segs.insert( bot.segs.begin(), D2<SBasis>(SBasis(Linear( f[X].at0()-2*tol*signx, f[X].at0())),
												  SBasis(Linear( f[Y].at0() )) ));
	top += Point(-tol*signx,  tol);
	bot += Point( tol*signx, -tol);

	if ( signy < 0 ){
		std::swap( top, bot );
		top += Point( 0,  2*tol);
		bot += Point( 0, -2*tol);
	}
	topside = make_cuts_independent(top);
	botside = make_cuts_independent(bot);
}


/*
 * Compute top and bottom boundaries of the L^infty nbhd of the graph of a *monotonic* function f.
 * if f is increasing, it is given by [f(t-tol)-tol, f(t+tol)+tol].
 * if not, it is [f(t+tol)-tol, f(t-tol)+tol].
 */
static
void computeLinfinityNeighborhood( Piecewise<SBasis> const &f, double tol, Piecewise<SBasis> &top, Piecewise<SBasis> &bot){
	top = f + tol;
	top.offsetDomain( - tol );
	top.cuts.insert( top.cuts.end(), f.domain().max() + tol);
	top.segs.insert( top.segs.end(), SBasis(Linear( f.lastValue() + tol )) );

	bot = f - tol;
	bot.offsetDomain( tol );
	bot.cuts.insert( bot.cuts.begin(), f.domain().min() - tol);
	bot.segs.insert( bot.segs.begin(), SBasis(Linear( f.firstValue() - tol )) );

    if (f.firstValue() > f.lastValue()) {
        std::swap(top, bot);
        top += 2 * tol;
        bot -= 2 * tol;
    }
}

/*
 * Returns the intervals over which the curve keeps its slope
 * in one of the 8 sectors delimited by x=0, y=0, y=x, y=-x.
 */
std::vector<Interval> monotonicSplit(D2<SBasis> const &p){
	std::vector<Interval> result;

	D2<SBasis> v = derivative(p);

	std::vector<double> someroots;
	std::vector<double> cuts (2,0.);
	cuts[1] = 1.;

	someroots = roots(v[X]);
	cuts.insert( cuts.end(), someroots.begin(), someroots.end() );

	someroots = roots(v[Y]);
	cuts.insert( cuts.end(), someroots.begin(), someroots.end() );

	//we could split in the middle to avoid computing roots again...
	someroots = roots(v[X]-v[Y]);
	cuts.insert( cuts.end(), someroots.begin(), someroots.end() );

	someroots = roots(v[X]+v[Y]);
	cuts.insert( cuts.end(), someroots.begin(), someroots.end() );

	sort(cuts.begin(),cuts.end());
	unique(cuts.begin(), cuts.end() );

	for (unsigned i=1; i<cuts.size(); i++){
		result.push_back( Interval( cuts[i-1], cuts[i] ) );
	}
	return result;
}

//std::vector<Interval> level_set( D2<SBasis> const &f, Rect region){
//	std::vector<Interval> x_in_reg = level_set( f[X], region[X] );
//	std::vector<Interval> y_in_reg = level_set( f[Y], region[Y] );
//	std::vector<Interval> result = intersect ( x_in_reg, y_in_reg );
//	return result;
//}

/*TODO: remove this!!!
 * the minimum would be to move it to piecewise.h but this would be stupid.
 * The best would be to let 'compose' be aware of extension modes (constant, linear, polynomial..)
 * (I think the extension modes (at start and end) should be properties of the pwsb).
 */
static
void prolongateByConstants( Piecewise<SBasis> &f, double paddle_width ){
	if ( f.size() == 0 ) return; //do we have a covention about the domain of empty pwsb?
	f.cuts.insert( f.cuts.begin(), f.cuts.front() - paddle_width );
	f.segs.insert( f.segs.begin(), SBasis( f.segs.front().at0() ) );
	f.cuts.insert( f.cuts.end(), f.cuts.back() + paddle_width );
	f.segs.insert( f.segs.end(), SBasis( f.segs.back().at1() ) );
}

static
bool compareIntersectionsTimesX( SmashIntersection const &inter1, SmashIntersection const &inter2 ){
	return inter1.times[X].min() < inter2.times[Y].min();
}
/*Fuse contiguous intersection domains
 *
 */
static
void cleanup_and_fuse( std::vector<SmashIntersection> &inters ){
	std::sort( inters.begin(), inters.end(), compareIntersectionsTimesX);
	for (unsigned i=0; i < inters.size(); i++ ){
		for (unsigned j=i+1; j < inters.size() && inters[i].times[X].intersects( inters[j].times[X]) ; j++ ){
			if (inters[i].times[Y].intersects( inters[j].times[Y] ) ){
				inters[i].times.unionWith(inters[j].times);
				inters[i].bbox.unionWith(inters[j].bbox);
				inters.erase( inters.begin() + j );
			}
		}
	}
}

/* Computes the intersection of two sets given as (ordered) union intervals.
 */
static
std::vector<Interval> intersect( std::vector<Interval> const &a, std::vector<Interval> const &b){
	std::vector<Interval> result;
	//TODO: use order to optimize this!
	for (auto i : a){
		for (auto j : b){
			OptInterval c( i );
			c &= j;
			if ( c ) {
				result.push_back( *c );
			}
		}
	}
	return result;
}

/* Returns the intervals over which the curves are in the
 * tol-neighborhood one of the other for the L_\infty metric.
 * WARNING: each curve is supposed to be a graph over x or y axis
 * (but not necessarily the same axis for both) and the smaller
 * the slope the better (typically <=45°).
 */
std::vector<SmashIntersection> monotonic_smash_intersect( D2<SBasis> const &a, D2<SBasis> const &b, double tol){
	using std::swap;

	// a and b or X and Y may have to be exchanged, so make local copies.
	D2<SBasis> aa = a;
	D2<SBasis> bb = b;
	bool swapresult = false;
	bool swapcoord = false;//debug only!

	//if the (enlarged) bounding boxes don't intersect, stop.
	OptRect abounds = bounds_fast( a );
	OptRect bbounds = bounds_fast( b );
	if ( !abounds || !bbounds ) return std::vector<SmashIntersection>();
	abounds->expandBy(tol);
	if ( !(abounds->intersects(*bbounds))){
		return std::vector<SmashIntersection>();
	}

	//Choose the best curve to be re-parametrized by x or y values.
	OptRect dabounds = bounds_exact(derivative(a));
	OptRect dbbounds = bounds_exact(derivative(b));
	if	( dbbounds->min().length() > dabounds->min().length() ){
		aa=b;
		bb=a;
		swap( dabounds, dbbounds );
		swapresult = true;
	}

	//Choose the best coordinate to use as new parameter
	double dxmin = std::min( std::abs((*dabounds)[X].max()), std::abs((*dabounds)[X].min()) );
	double dymin = std::min( std::abs((*dabounds)[Y].max()), std::abs((*dabounds)[Y].min()) );
	if ( (*dabounds)[X].max()*(*dabounds)[X].min() < 0 ) dxmin=0;
	if ( (*dabounds)[Y].max()*(*dabounds)[Y].min() < 0 ) dymin=0;
	assert (dxmin>=0 && dymin>=0);

	if (dxmin < dymin) {
		aa = D2<SBasis>( aa[Y], aa[X] );
		bb = D2<SBasis>( bb[Y], bb[X] );
		swapcoord = true;
	}

	//re-parametrize aa by the value of x.
	Interval x_range_strict( aa[X].at0(), aa[X].at1() );
	Piecewise<SBasis> y_of_x = pw_compose_inverse(aa[Y],aa[X], 2, 1e-5);

	//Compute top and bottom boundaries of the L^infty nbhd of aa.
	Piecewise<SBasis> top_ay, bot_ay;
	computeLinfinityNeighborhood( y_of_x, tol, top_ay, bot_ay);

	Interval ax_range = top_ay.domain();//i.e. aa[X] domain ewpanded by tol.
	std::vector<Interval> bx_in_ax_range = level_set(bb[X], ax_range );

	// find times when bb is in the neighborhood of aa.
	std::vector<Interval> tbs;
	for (auto & i : bx_in_ax_range){
		D2<Piecewise<SBasis> > bb_in;
		bb_in[X] = Piecewise<SBasis> ( portion( bb[X], i ) );
		bb_in[Y] = Piecewise<SBasis> ( portion( bb[Y], i) );
		bb_in[X].setDomain( i );
		bb_in[Y].setDomain( i );

		Piecewise<SBasis> h;
		Interval level;
		h = bb_in[Y] - compose( top_ay, bb_in[X] );
		level = Interval( -infinity(), 0 );
		std::vector<Interval> rts_lo = level_set( h, level);
		h = bb_in[Y] - compose( bot_ay, bb_in[X] );
		level = Interval( 0, infinity());
		std::vector<Interval> rts_hi = level_set( h, level);

		std::vector<Interval> rts = intersect( rts_lo, rts_hi );
		tbs.insert(tbs.end(), rts.begin(),  rts.end()  );
	}

	std::vector<SmashIntersection> result(tbs.size(), SmashIntersection());

	/* for each solution I, find times when aa is in the neighborhood of bb(I).
	 * (Note: the preimage of bb[X](I) by aa[X], enlarged by tol, is a good approximation of this:
	 * it would give points in the 2*tol neighborhood of bb (if the slope of aa is never more than 1).
	 *  + faster computation.
	 *  - implies little jumps depending on the subdivision of the input curve into monotonic pieces
	 *  and on the choice of preferred axis. If noticeable, these jumps would feel random to the user :-(
	 */
	for (unsigned j=0; j<tbs.size(); j++){
		result[j].times[Y] = tbs[j];
		std::vector<Interval> tas;
		//TODO: replace this by some option in the "compose(pw,pw)" method!
		Piecewise<SBasis> fat_y_of_x = y_of_x;
		prolongateByConstants( fat_y_of_x, 100*(1+tol) );

		D2<Piecewise<SBasis> > top_b, bot_b;
		D2<SBasis> bbj = portion( bb, tbs[j] );
		computeLinfinityNeighborhood( bbj, tol, top_b, bot_b );

		Piecewise<SBasis> h;
		Interval level;
		h = top_b[Y] - compose( fat_y_of_x, top_b[X] );
		level = Interval( +infinity(), 0 );
		std::vector<Interval> rts_top = level_set( h, level);
		for (auto & idx : rts_top){
			idx = Interval( top_b[X].valueAt( idx.min() ),
									 top_b[X].valueAt( idx.max() ) );
		}
		assert( rts_top.size() == 1 );

		h = bot_b[Y] - compose( fat_y_of_x, bot_b[X] );
		level = Interval( 0, -infinity());
		std::vector<Interval> rts_bot = level_set( h, level);
		for (auto & idx : rts_bot){
			idx = Interval( bot_b[X].valueAt( idx.min() ),
									 bot_b[X].valueAt( idx.max() ) );
		}
		assert( rts_bot.size() == 1 );
		rts_top = intersect( rts_top, rts_bot );
		assert (rts_top.size() == 1);
		Interval x_dom = rts_top[0];

		if ( x_dom.max() <= x_range_strict.min() ){
			tas.push_back( Interval ( ( aa[X].at0() < aa[X].at1() ) ? 0 : 1 ) );
		}else if ( x_dom.min() >= x_range_strict.max() ){
			tas.push_back( Interval ( ( aa[X].at0() < aa[X].at1() ) ? 1 : 0 ) );
		}else{
			tas = level_set(aa[X], x_dom );
		}
		assert( tas.size()==1 );
		result[j].times[X] = tas.front();

		result[j].bbox = Rect( bbj.at0(), bbj.at1() );
		Interval y_dom( aa[Y](result[j].times[X].min()), aa[Y](result[j].times[X].max()) );
		result[j].bbox.unionWith( Rect( x_dom, y_dom ) );
	}

	if (swapresult) {
		for (auto & i : result){
			swap( i.times[X], i.times[Y]);
		}
	}
	if (swapcoord) {
		for (auto & i : result){
			swap( i.bbox[X], i.bbox[Y] );
		}
	}

	//TODO: cleanup result? fuse contiguous intersections?
	return result;
}

std::vector<SmashIntersection> smash_intersect( D2<SBasis> const &a, D2<SBasis> const &b, double tol){
	std::vector<SmashIntersection> result;

	std::vector<Interval> acuts = monotonicSplit(a);
	std::vector<Interval> bcuts = monotonicSplit(b);
	for (auto & acut : acuts){
		D2<SBasis> ai = portion( a, acut);
		for (auto & bcut : bcuts){
			D2<SBasis> bj = portion( b, bcut);
			std::vector<SmashIntersection> ai_cap_bj = monotonic_smash_intersect( ai, bj, tol );
			for (auto & k : ai_cap_bj){
				k.times[X] = k.times[X] * acut.extent() + acut.min();
				k.times[Y] = k.times[Y] * bcut.extent() + bcut.min();
			}
			result.insert( result.end(), ai_cap_bj.begin(), ai_cap_bj.end() );
		}
	}
	cleanup_and_fuse( result );
	return result;
}

}

/*
  Local Variables:
  mode:c++
  c-file-style:"stroustrup"
  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
  indent-tabs-mode:nil
  fill-column:99
  End:
*/
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :