//============================================================================
//
// This file is part of the 'reorientmeshfaces' project.
//
// This software is covered by the following BSD license, except for portions
// derived from other works which are covered by their respective licenses.
// For full licensing information including reproduction of these external
// licenses, see the file LICENSE.txt provided in the documentation.
//
// Copyright (c) 2015, Evangelos Kalogerakis, UMass Amherst
//
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
//
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// * Neither the name of the copyright holders nor the names of contributors
// to this software may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
//============================================================================

#include "Mesh.h"


// constructor for cleaning up / preparing an input mesh for shape sampling
Mesh::Mesh(const string& filename, bool add_faces, bool do_nothing)
{
	mesh = TriMesh::read(filename.c_str());  
	if (mesh == NULL)
		return;
	if (mesh->faces.size() == 0) 
	{
		mesh = NULL;
		std::cout << "Mesh has 0 faces! Skipping this mesh..." << std::endl;
		return;
	}
  if (do_nothing)
    return;

	tree = NULL;
	mesh->bbox.valid = false; mesh->need_bsphere();
	
	for (size_t i = 0; i < mesh->vertices.size(); i++) 
	{
		mesh->vertices[i] /= mesh->bsphere.r;
	}


	mesh->normals.clear();
	mesh->need_normals();		// this compute face areas as well!
	mesh->faces.reserve( 2 * mesh->faces.size() );

	for (size_t i = 0; i < mesh->faces.size(); i++) 
	{
		Vector3 v[3] = { Vector3( mesh->vertices[mesh->faces[i][0]][0], mesh->vertices[mesh->faces[i][0]][1], mesh->vertices[mesh->faces[i][0]][2] ),
			Vector3( mesh->vertices[mesh->faces[i][1]][0], mesh->vertices[mesh->faces[i][1]][1], mesh->vertices[mesh->faces[i][1]][2] ),
			Vector3( mesh->vertices[mesh->faces[i][2]][0], mesh->vertices[mesh->faces[i][2]][1], mesh->vertices[mesh->faces[i][2]][2] ) };
		MyCustomTriangle tri(MyCustomTriangleVertexTriple(v[0], v[1], v[2], i));
		triangles.push_back(tri);
	}
	tree = new TKDTree(triangles.begin(),triangles.end());
	


	size_t num_original_faces = mesh->faces.size();
	for (size_t i = 0; i < num_original_faces; i++) 
	{
		if (i % 100 == 0)
			std::cout << "Processed " << i+1 << " out of " << num_original_faces << std::endl;
		
		vec v = mesh->faces[i].facenormal;
		vec u = mesh->vertices[mesh->faces[i][1]] - mesh->vertices[mesh->faces[i][0]];
		u = normalize(u);
		u = _normalize(u CROSS v);

		int hitsP = 0;
		int numPrays = 0;
		vec tv, ttv;
		float angle;
		bool hit;		
#pragma omp parallel for private(angle, tv, ttv, hit)
		for (int a = 0; a <= 3; a++) 
		{ 
			angle = float(a) * 20.0f;
			tv = _normalize( cosd(angle)*v + (1 - cosd(angle))*(v DOT u)*u + sind(angle)*(u CROSS v) );
			for (float rotangle = 0; rotangle < 360; rotangle += 45.0f) 
			{
				ttv = _normalize( cosd(rotangle)*tv + (1 - cosd(rotangle))*(tv DOT v)*v + sind(rotangle)*(v CROSS tv) );
				hit = rayShoot(ttv, mesh->faces[i].faceCenter, -1 ) > 0.0f;
#pragma omp atomic
				hitsP += hit;
#pragma omp atomic
				numPrays++;
			}
		}

		int hitsN = 0;
		int numNrays = 0;
		v = -mesh->faces[i].facenormal;

#pragma omp parallel for private(angle, tv, ttv, hit)
		for (int a = 0; a <= 3; a++) 
		{
			angle = float(a) * 20.0f;
			tv = _normalize( cosd(angle)*v + (1 - cosd(angle))*(v DOT u)*u + sind(angle)*(u CROSS v) );
			for (float rotangle = 0; rotangle < 360; rotangle += 45.0f) {				
				ttv = _normalize( cosd(rotangle)*tv + (1 - cosd(rotangle))*(tv DOT v)*v + sind(rotangle)*(v CROSS tv) );
				hit = rayShoot(ttv, mesh->faces[i].faceCenter, -1) > 0.0f;
#pragma omp atomic
				hitsN += hit;
#pragma omp atomic
				numNrays++;
				//hitsN += (rayShoot(ttv, mesh->vertices[mesh->faces[i][0]], -1) > 0.0f);
				//hitsN += (rayShoot(ttv, mesh->vertices[mesh->faces[i][1]], -1) > 0.0f);
				//hitsN += (rayShoot(ttv, mesh->vertices[mesh->faces[i][2]], -1) > 0.0f);
			}
		}
		assert( numNrays == numPrays );



		if (hitsN < hitsP) 
		{
			//std::cout << "Inverting face normals...";
			std::swap( mesh->faces[i][0], mesh->faces[i][2] );
			mesh->faces[i].facenormal = -mesh->faces[i].facenormal;
		} 
		float max_ratio_hit = max( (float)hitsN / (float)numNrays , (float)hitsP / (float)numPrays);
		float min_ratio_hit = min( (float)hitsN / (float)numNrays , (float)hitsP / (float)numPrays);

    if (add_faces)
    {
      if (min_ratio_hit > .9f || max_ratio_hit < .95f) // if the "least visible" orientation has less than 95% of the rays intersecting, then create a double face
      {
        int num_verts = mesh->vertices.size();
        mesh->vertices.push_back(mesh->vertices[mesh->faces[i][0]] - (1e-6f) * mesh->faces[i].facenormal);
        mesh->vertices.push_back(mesh->vertices[mesh->faces[i][1]] - (1e-6f) * mesh->faces[i].facenormal);
        mesh->vertices.push_back(mesh->vertices[mesh->faces[i][2]] - (1e-6f) * mesh->faces[i].facenormal);
        mesh->faces.push_back(TriMesh::Face(num_verts + 2, num_verts + 1, num_verts));
      }
    }


		if ( i == num_original_faces-1)
			std::cout << "Processed " << i+1 << " out of " << num_original_faces << ", added " << mesh->faces.size() - num_original_faces << " faces." << std::endl;

	}
	triangles.clear();
	delete tree;
	tree = NULL;

	std::cout << " Done." << std::endl;
	std::cout << std::endl;






	/////////////////////////////////////////////////////////
	//// subdivide
	/////////////////////////////////////////////////////////
	//while (mesh->faces.size() < numsamples )
	//	subdivSelectively(mesh, 0.0f);

	//std::cout << std::endl;

	//for (int subdiv_iter = 1; subdiv_iter<=10; subdiv_iter++)
	//{
	//	vector<float> maximum_length;
	//	mesh->normals.clear();
	//	mesh->need_normals();
	//	int nf = mesh->faces.size();
	//	for (int i = 0; i < nf; i++)
	//	{
	//		vec e1 = mesh->vertices[mesh->faces[i][0]] - mesh->vertices[mesh->faces[i][1]];
	//		vec e2 = mesh->vertices[mesh->faces[i][0]] - mesh->vertices[mesh->faces[i][2]];
	//		vec e3 = mesh->vertices[mesh->faces[i][1]] - mesh->vertices[mesh->faces[i][2]];
	//		maximum_length.push_back( max( max( len(e1), len(e2) ), len(e3) ) );
	//	}
	//	float max_edge_length_thr = 3.333f * median( maximum_length );
	//	maximum_length.clear();

	//	std::cout << "Subdivision round " << subdiv_iter << ", #faces= " << mesh->faces.size() << " ...";
	//	if ( subdivSelectively(mesh, max_edge_length_thr) == false )
	//		break;
	//}

	//std::cout << std::endl;
	/////////////////////////////////////////////////////////



	/////////////////////////////////////////////////////////
	//// dereference redundant vertices
	/////////////////////////////////////////////////////////
	//if ( !VERY_LARGE_MESHES_MODE )
	//{
	//	std::cout << " Checking for duplicate vertices and welding when necessary... " << std::endl;

	//	int deref_round = 0;
	//	mesh->normals.clear();
	//	float feature_size = .01f * mesh->feature_size();
	//	while(deref_round <= 12)
	//	{
	//		deref_round++;
	//		bool deref = false;
	//		std::cout << " Round " << deref_round << ", ";

	//		vector<bool> destroyed_vertices( mesh->vertices.size(), false );
	//		vector<NamedPoint> points;
	//		for (size_t i = 0; i < mesh->vertices.size(); i++)
	//			points.push_back( NamedPoint(  mesh->vertices[i][0], mesh->vertices[i][1], mesh->vertices[i][2], i ) );
	//		PKDTree stree(points.begin(), points.end());
	//		for (size_t i = 0; i < mesh->vertices.size(); i++) 
	//		{
	//			if ( destroyed_vertices[ i ] )
	//				continue;
	//			size_t nearest_point_id;
	//			NamedPoint query( mesh->vertices[i][0], mesh->vertices[i][1], mesh->vertices[i][2] );
	//			BoundedSortedArray<PKDTree::Neighbor> k_closest_elems(2);
	//			stree.kClosestElements<MetricL2>(query, k_closest_elems, -1);  // -1 means there's no limit on the maximum allowed
	//			for (int kn = 0; kn < k_closest_elems.size(); kn++)
	//			{
	//				if ( stree.getElements()[k_closest_elems[kn].getIndex()].id != i) 
	//				{
	//					nearest_point_id = stree.getElements()[k_closest_elems[kn].getIndex()].id;
	//					break;
	//				}
	//			}
	//			if ( dist2(mesh->vertices[i], mesh->vertices[nearest_point_id]) < sqr(feature_size) ) 
	//			{
	//				deref = true;
	//				destroyed_vertices[ nearest_point_id ] = true;
	//				for (size_t f = 0; f < mesh->faces.size(); f++) 
	//				{
	//					if (mesh->faces[f][0] == nearest_point_id)
	//						mesh->faces[f][0] = i;
	//					if (mesh->faces[f][1] == nearest_point_id)
	//						mesh->faces[f][1] = i;
	//					if (mesh->faces[f][2] == nearest_point_id)
	//						mesh->faces[f][2] = i;              
	//				}
	//			}
	//		}
	//		destroyed_vertices.clear();
	//		points.clear();
	//		remove_unused_vertices(mesh);
	//		if (deref == false)
	//			break;
	//	}
	//	remove_sliver_faces(mesh);
	//	std::cout << std::endl;
	//}
} // end of correct_mesh_mode


void Mesh::write(const string& filename)
{
	std::cout << "Writing " << filename << std::endl;
	mesh->write(filename.c_str(), true);
}


float Mesh::rayShoot(vec& line, point& origin, float limdist) 
{
	float eps = 1e-5f;
	Ray3 ray( Vector3(origin[0] + eps * line[0], origin[1] + eps * line[1], origin[2] + eps * line[2]), Vector3(line[0], line[1], line[2]) );
	//Real hit_time = tree->rayIntersectionTime(ray, limdist + eps);
	Real hit_time = tree->rayIntersectionTime(ray);
	if (hit_time >= 0)
		return hit_time;
	else
		return -1.0f;

	//list<Primitive_id> primitives;
	//tree->all_intersected_primitives(ray_query, std::back_inserter(primitives));
	//float currentDist = FLT_MAX;


	//for(list<Primitive_id>::iterator i = primitives.begin(); i != primitives.end(); ++i) {
	//  Primitive_id primitive_id = *i; 
	//  size_t t = primitive_id->id;    

	//  if (t == f)
	//    continue;

	//  if ( (line DOT mesh->faces[t].facenormal) < 0 )
	//    continue;

	//  vec e1 = mesh->vertices[mesh->faces[t][1]] - mesh->vertices[mesh->faces[t][0]];
	//  vec e2 = mesh->vertices[mesh->faces[t][2]] - mesh->vertices[mesh->faces[t][0]];
	//  vec p = line CROSS e2;
	//  float tmp = p DOT e1;

	//  if ( tmp < EPSILON && tmp > -EPSILON) 
	//    continue;

	//  tmp = 1.0f / tmp;
	//  vec s = origin - mesh->vertices[mesh->faces[t][0]];
	//  float u = tmp * (s DOT p);
	//  if (u < 0.0 || u > 1.0)
	//    continue;

	//  vec q = s CROSS e1;
	//  float v = tmp * (line DOT q);
	//  if (v < 0.0 || v > 1.0 || u + v > 1.0)
	//    continue;


	//  float distance = tmp * (e2 DOT q);		
	//  if (distance <= 0.0f) 
	//    continue;

	//  if (distance < currentDist)
	//    currentDist = distance;
	//} // end of searching of intersections

	//if ( currentDist == FLT_MAX ) {
	//  currentDist = -1.0f;
	//}
	//primitives.clear();

	//return currentDist;
}




