// This file is part of libigl, a simple c++ geometry processing library.
// 
// Copyright (C) 2016 Yotam Gingold <yotam@yotamgingold.com>
// 
// 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 "seam_edges.h"
#include <unordered_map>
#include <unordered_set>
#include <cassert>

namespace {
    // Computes the orientation of `c` relative to the line between `a` and `b`.
    // Assumes 2D vector input.
    // Based on: https://www.cs.cmu.edu/~quake/robust.html
    template< typename scalar_t >
    inline scalar_t orientation(
        const Eigen::Matrix< scalar_t, 2, 1 >& a,
        const Eigen::Matrix< scalar_t, 2, 1 >& b,
        const Eigen::Matrix< scalar_t, 2, 1 >& c
        ) {
        typedef Eigen::Matrix< scalar_t, 2, 1 > Vector2S;
        const Vector2S row0 = a - c;
        const Vector2S row1 = b - c;
        return row0(0)*row1(1) - row1(0)*row0(1);
    }
}

// I have verified that this function produces the exact same output as
// `find_seam_fast.py` for `cow_triangled.obj`.
template <typename DerivedV, typename DerivedF, typename DerivedT>
IGL_INLINE void seam_edges(
    const Eigen::PlainObjectBase<DerivedV>& V,
    const Eigen::PlainObjectBase<DerivedT>& TC,
    const Eigen::PlainObjectBase<DerivedF>& F,
    const Eigen::PlainObjectBase<DerivedF>& FTC,
    Eigen::PlainObjectBase<DerivedF>& seams,
    Eigen::PlainObjectBase<DerivedF>& boundaries,
    Eigen::PlainObjectBase<DerivedF>& foldovers
    )
{
    // Assume triangles.
    assert( F.cols() == 3 );
    assert( F.cols() == FTC.cols() );
    assert( F.rows() == FTC.rows() );
    
    // Assume 2D texture coordinates (foldovers tests).
    assert( TC.cols() == 2 );
    typedef Eigen::Matrix< typename DerivedT::Scalar, 2, 1 > Vector2S;
    
    seams     .setZero( 3*F.rows(), 4 );
    boundaries.setZero( 3*F.rows(), 2 );
    foldovers .setZero( 3*F.rows(), 4 );
    
    int num_seams = 0;
    int num_boundaries = 0;
    int num_foldovers = 0;
    
    // A map from a pair of vertex indices to the index (face and endpoints)
    // into face_position_indices.
    // The following should be true for every key, value pair:
    //    key == face_position_indices[ value ]
    // This gives us a "reverse map" so that we can look up other face
    // attributes based on position edges.
    // The value are written in the format returned by numpy.where(),
    // which stores multi-dimensional indices such as array[a0,b0], array[a1,b1]
    // as ( (a0,a1), (b0,b1) ).
    
    // We need to make a hash function for our directed edges.
    // We'll use i*V.rows() + j.
    typedef std::pair< typename DerivedF::Scalar, typename DerivedF::Scalar > directed_edge;
	const int numV = V.rows();
	const int numF = F.rows();
	auto edge_hasher = [numV]( directed_edge const& e ) { return e.first*numV + e.second; };
	// When we pass a hash function object, we also need to specify the number of buckets.
	// The Euler characteristic says that the number of undirected edges is numV + numF -2*genus.
	std::unordered_map< directed_edge, std::pair< int, int >, decltype( edge_hasher ) > directed_position_edge2face_position_index( 2*( numV + numF ), edge_hasher );
    for( int fi = 0; fi < F.rows(); ++fi ) {
        for( int i = 0; i < 3; ++i ) {
            const int j = ( i+1 ) % 3;
            directed_position_edge2face_position_index[ std::make_pair( F(fi,i), F(fi,j) ) ] = std::make_pair( fi, i );
        }
    }
    
    // First find all undirected position edges (collect both orientations of
    // the directed edges).
    std::unordered_set< directed_edge, decltype( edge_hasher ) > undirected_position_edges( numV + numF, edge_hasher );
    for( const auto& el : directed_position_edge2face_position_index ) {
        undirected_position_edges.insert( el.first );
        undirected_position_edges.insert( std::make_pair( el.first.second, el.first.first ) );
    }
    
    // Now we will iterate over all position edges.
    // Seam edges are the edges whose two opposite directed edges have different
    // texcoord indices (or one doesn't exist at all in the case of a mesh
    // boundary).
    for( const auto& vp_edge : undirected_position_edges ) {
        const auto vp_edge_reverse = std::make_pair( vp_edge.second, vp_edge.first );
        // If it and its opposite exist as directed edges, check if their
        // texture coordinate indices match.
        if( directed_position_edge2face_position_index.count( vp_edge ) &&
            directed_position_edge2face_position_index.count( vp_edge_reverse ) ) {
            const auto forwards = directed_position_edge2face_position_index[ vp_edge ];
            const auto backwards = directed_position_edge2face_position_index[ vp_edge_reverse ];
            
            // NOTE: They should never be equal.
            assert( forwards != backwards );
            
            // NOTE: Non-matching seam edges will show up twice, once as
            //       ( forwards, backwards ) and once as
            //       ( backwards, forwards ). We don't need to examine both,
            //       so continue only if forwards < backwards.
            if( forwards < backwards ) continue;

            // If the texcoord indices match (are similarly flipped),
            // this edge is not a seam. It could be a foldover.
            if( std::make_pair( FTC( forwards.first, forwards.second ), FTC( forwards.first, ( forwards.second+1 ) % 3 ) )
                ==
                std::make_pair( FTC( backwards.first, ( backwards.second+1 ) % 3 ), FTC( backwards.first, backwards.second ) )
                ) {
                
                // Check for foldovers in UV space.
                // Get the edge (a,b) and the two opposite vertices's texture coordinates.
                const Vector2S a = TC.row( FTC( forwards.first,  forwards.second ) );
                const Vector2S b = TC.row( FTC( forwards.first, (forwards.second+1) % 3 ) );
                const Vector2S c_forwards  = TC.row( FTC( forwards .first, (forwards .second+2) % 3 ) );
                const Vector2S c_backwards = TC.row( FTC( backwards.first, (backwards.second+2) % 3 ) );
                // If the opposite vertices' texture coordinates fall on the same side
                // of the edge, we have a UV-space foldover.
                const auto orientation_forwards = orientation( a, b, c_forwards );
                const auto orientation_backwards = orientation( a, b, c_backwards );
                if( ( orientation_forwards > 0 && orientation_backwards > 0 ) ||
                    ( orientation_forwards < 0 && orientation_backwards < 0 )
                    ) {
                    foldovers( num_foldovers, 0 ) = forwards.first;
                    foldovers( num_foldovers, 1 ) = forwards.second;
                    foldovers( num_foldovers, 2 ) = backwards.first;
                    foldovers( num_foldovers, 3 ) = backwards.second;
                    num_foldovers += 1;
                }
            }
            
            // Otherwise, we have a non-matching seam edge.
            else {
                seams( num_seams, 0 ) = forwards.first;
                seams( num_seams, 1 ) = forwards.second;
                seams( num_seams, 2 ) = backwards.first;
                seams( num_seams, 3 ) = backwards.second;
                num_seams += 1;
            }
        }
        // Otherwise, the edge and its opposite aren't both in the directed
        // edges. One of them should be.
        else if( directed_position_edge2face_position_index.count( vp_edge ) ) {
            const auto forwards = directed_position_edge2face_position_index[ vp_edge ];
            boundaries( num_boundaries, 0 ) = forwards.first;
            boundaries( num_boundaries, 1 ) = forwards.second;
            num_boundaries += 1;
        }
        else if( directed_position_edge2face_position_index.count( vp_edge_reverse ) ) {
            const auto backwards = directed_position_edge2face_position_index[ vp_edge_reverse ];
            boundaries( num_boundaries, 0 ) = backwards.first;
            boundaries( num_boundaries, 1 ) = backwards.second;
            num_boundaries += 1;
        }
        else {
            // This should never happen! One of these two must have been seen.
            assert(
                directed_position_edge2face_position_index.count( vp_edge ) ||
                directed_position_edge2face_position_index.count( vp_edge_reverse )
                );
        }
    }
    
    seams     .conservativeResize( num_seams,      Eigen::NoChange_t() );
    boundaries.conservativeResize( num_boundaries, Eigen::NoChange_t() );
    foldovers .conservativeResize( num_foldovers,  Eigen::NoChange_t() );
}