/* This file is part of the source code for the following publication:
 * Assembling Self-Supporting Structures, Deuss et al., SIGGRAPH Asia 2014
 *
 * Copyright (C) 2014 Mario Deuss <mario.deuss@epfl.ch>
 *
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */

#include "Utils.h"

namespace SelfAssembly {

bool Utils::intersects(const V &a, const V &b, const V &c, 
						  const V &o, const V &dir, V* bc/*=NULL*/, 
						  double tMin /*=10e-6*/, double tMax /*= std::numeric_limits<double>::max()*/)
{
	//Barycentric coordinates
	V bcc=V::Zero();
	
	//taken from pharr&humphreys physically based rendering
	V  e1 = b - a;
	V  e2 = c - a;
	V  s1 = dir.cross(e2);
	
	double div = s1.dot(e1);
	if (div==0) return false; // no intersection due to degenerate triangle
	
	double inv = 1.0/div;
	
	// compute first barycentric coordinate
	V  dist = o - a;
	bcc(1) = dist.dot(s1) * inv;
	
	
	// compute second barycentric coordinate
	V  s2 = dist.cross(e1);
	bcc(2) = dir.dot(s2) * inv;
	
	//check if below max t
	double tInt = e2.dot(s2) *  inv;
	if(tInt > tMax || tInt < tMin){
		return false;
	}
	
	// compute ray plane intersection
	//V pInt = a*(1-bc(1)-bc(2)) + b*bc(1) + c*bc(2);
	bcc(0)=1-bcc(1)-bcc(2);
	
	if ((bcc(1)<0.0) || (bcc(1)>1.0)) {
		return false; // no intersection
	}
	else if ((bcc(2)<0.0) || (bcc(1)+bcc(2)>1.0)) {
		return false; // no intersection
	}
	else {
		if(bc!=NULL){(*bc)=bcc;}
		return true;
	}
}

V Utils::sampleCircle()
{
	double theta = 2. * M_PI *(double)rand()/(double)RAND_MAX;
	double r = std::sqrt( (double)rand()/(double)RAND_MAX);
	V v; v<< r*std::cos(theta) , r*std::sin(theta), 0.;
	return v;
}

V Utils::incircleCenter(const V &a, const V &b, const V &c)
{
	double la = (b-c).norm();
	double lb = (c-a).norm();
	double lc = (a-b).norm();
	return (la*a + lb*b + lc*c)/(la+lb+lc);
}

double Utils::incircleRadius(const V &a, const V &b, const V &c)
{
	double la = (b-c).norm();
	double lb = (c-a).norm();
	double lc = (a-b).norm();
	return 0.5* std::sqrt(((lb+lc-la)*(lc+la-lb)*(la+lb-lc))/(la+lb+lc));
}


}
