//==============================================================================
/*! \file
 * OpenMesh Toolkit for mesh analysis    \n
 * Copyright (c) 2010 by Rostislav Hulik     \n
 *
 * Author:  Rostislav Hulik, rosta.hulik@gmail.com  \n
 * Date:    2010/11/21                          \n
 *
 * This file is part of software developed for support of Rostislav Hulik's dissertation thesis at dcgm-robotics@FIT group.
 *
 * This file is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This file is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this file.  If not, see <http://www.gnu.org/licenses/>.
 * 
 * Description:
 * - Class smoothes the mesh using the Local Projections method
 */
 
#ifndef _OMSMOOTHER_H_
#define _OMSMOOTHER_H_

#include <OMToolkit\OMTypes.h>
#include <OMToolkit\OMGaussFunc.h>
#include <OpenMesh\Tools\Smoother\JacobiLaplaceSmootherT.hh>
namespace OMToolkit {

class OMSmoother
{
	public:
		typedef Types::ModuleMeshd MeshT;
		typedef MeshT::Point PointT;
		typedef MeshT::Normal NormalT;
		typedef MeshT::Scalar ScalarT;
		typedef MeshT::VertexHandle VertexH;

		OMSmoother(MeshT *mesh, OpenMesh::VPropHandleT<NormalT> &curvVectorHandle, OpenMesh::VPropHandleT<ScalarT> &curvHandle, ScalarT gaussF);
		~OMSmoother();

		void smooth(int n, ScalarT alpha);

	private:

		void computeWeightsAround(VertexH vertex);
		OpenMesh::VPropHandleT<PointT>	original_positions;
		OpenMesh::VPropHandleT<NormalT> original_normals;
		OpenMesh::VPropHandleT<ScalarT>	weight;
		OpenMesh::VPropHandleT<NormalT> weight_direction;
		OMGaussFunction<ScalarT> *gaussFunc;
		MeshT *m_mesh;
};



typedef Types::ModuleMeshd MeshT;

class smoother2 : public OpenMesh::Smoother::JacobiLaplaceSmootherT<MeshT>
{
public:
	

	smoother2(MeshT *mesh) : OpenMesh::Smoother::JacobiLaplaceSmootherT<Types::ModuleMeshd>::JacobiLaplaceSmootherT(*mesh)
	{
		
	}
	typedef SmootherT<MeshT>                   Base;
	enum LaplaceWeighting {jedna};
	void compute_weights(LaplaceWeighting _weighting)
	{
		MeshT::VertexIter        v_it, v_end(Base::mesh_.vertices_end());
		MeshT::EdgeIter          e_it, e_end(Base::mesh_.edges_end());
		 MeshT::HalfedgeHandle    heh0, heh1, heh2;
		 MeshT::VertexHandle      v0, v1;
		//const  MeshT::Point       *p0, *p1, *p2;
		 MeshT::Normal            d0, d1;
		MeshT::Scalar            weight, lb(-1.0), ub(1.0);



		OpenMesh::VPropHandleT<float> vertex_weights_;
		Base::mesh_.add_property(vertex_weights_);

		OpenMesh::EPropHandleT<float> edge_weights_;
		Base::mesh_.add_property(edge_weights_);
		// init vertex weights
		for (v_it=Base::mesh_.vertices_begin(); v_it!=v_end; ++v_it)
			Base::mesh_.property(vertex_weights_, v_it) = 0.0;

		
    // Uniform weighting
    
      for (e_it=Base::mesh_.edges_begin(); e_it!=e_end; ++e_it)
      {
		heh0   = Base::mesh_.halfedge_handle(e_it.handle(), 0);
		heh1   = Base::mesh_.halfedge_handle(e_it.handle(), 1);
		v0     = Base::mesh_.to_vertex_handle(heh0);
		v1     = Base::mesh_.to_vertex_handle(heh1);
	
		Base::mesh_.property(edge_weights_, e_it) = 1.0;
		Base::mesh_.property(vertex_weights_, v0) += 1.0;
		Base::mesh_.property(vertex_weights_, v1) += 1.0;
      }

    


   // // Cotangent weighting
   // case CotWeighting:
   // {
   //   for (e_it=Base::mesh_.edges_begin(); e_it!=e_end; ++e_it)
   //   {
			//weight = 0.0;
	
			//heh0   = Base::mesh_.halfedge_handle(e_it.handle(), 0);
			//v0     = Base::mesh_.to_vertex_handle(heh0);
			//p0     = &Base::mesh_.point(v0);
	
			//heh1   = Base::mesh_.halfedge_handle(e_it.handle(), 1);
			//v1     = Base::mesh_.to_vertex_handle(heh1);
			//p1     = &Base::mesh_.point(v1);
	
			//heh2   = Base::mesh_.next_halfedge_handle(heh0);
			//p2     = &Base::mesh_.point(Base::mesh_.to_vertex_handle(heh2));
			//d0     = (*p0 - *p2); d0.normalize();
			//d1     = (*p1 - *p2); d1.normalize();
			//weight += 1.0 / tan(acos(std::max(lb, std::min(ub, dot(d0,d1) ))));
	
			//heh2   = Base::mesh_.next_halfedge_handle(heh1);
			//p2     = &Base::mesh_.point(Base::mesh_.to_vertex_handle(heh2));
			//d0     = (*p0 - *p2); d0.normalize();
			//d1     = (*p1 - *p2); d1.normalize();
			//weight += 1.0 / tan(acos(std::max(lb, std::min(ub, dot(d0,d1) ))));
	
			//Base::mesh_.property(edge_weights_, e_it) = weight;
			//Base::mesh_.property(vertex_weights_, v0)  += weight;
			//Base::mesh_.property(vertex_weights_, v1)  += weight;
   //   }
   //   break;
   // }
  

  
  // invert vertex weights:
  // before: sum of edge weights
  // after: one over sum of edge weights
  for (v_it=Base::mesh_.vertices_begin(); v_it!=v_end; ++v_it)
  {
    weight = Base::mesh_.property(vertex_weights_, v_it);
    if (weight)
      Base::mesh_.property(vertex_weights_, v_it) = 1.0 / weight;
  }
	}
};


} // namespace
#endif