yig/seam_edges.h

## seam_edges.h
// Copyright (C) 2016 Yotam Gingold <yotam@yotamgingold.com>
// GitHub Gist: https://gist.github.com/yig/bea834b6e5262490611cdb4470a2f34c

/*
Boost Software License - Version 1.0 - August 17th, 2003

Permission is hereby granted, free of charge, to any person or organization
obtaining a copy of the software and accompanying documentation covered by
this license (the "Software") to use, reproduce, display, distribute,
execute, and transmit the Software, and to prepare derivative works of the
Software, and to permit third-parties to whom the Software is furnished to
do so, all subject to the following:

The copyright notices in the Software and this entire statement, including
the above license grant, this restriction and the following disclaimer,
must be included in all copies of the Software, in whole or in part, and
all derivative works of the Software, unless such copies or derivative
works are solely in the form of machine-executable object code generated by
a source language processor.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
DEALINGS IN THE SOFTWARE.
*/

#ifndef SEAM_EDGES_H
#define SEAM_EDGES_H

#include <Eigen/Core>
#include <unordered_map>
#include <unordered_set>
#include <cassert>

/*
    Returns all UV-space boundaries of a mesh.
    Inputs:
        V: The positions of the input mesh.
        TC: The 2D texture coordinates of the input mesh (2 columns).
        F: A manifold-with-boundary triangle mesh, as vertex indices into `V` for the three vertices of each triangle.
        FTC: Indices into `TC` for the three vertices of each triangle.
    Outputs:
        seams: Edges where the forwards and backwards directions have different texture
               coordinates, as a #seams-by-4 matrix of indices.
               Each row is organized as [ forward_face_index, forward_face_vertex_index,
               backwards_face_index, backwards_face_vertex_index ]
               such that one side of the seam is the edge:
                    F[ seams( i, 0 ), seams( i, 1 ) ], F[ seams( i, 0 ), (seams( i, 1 ) + 1) % 3 ]
               and the other side is the edge:
                    F[ seams( i, 2 ), seams( i, 3 ) ], F[ seams( i, 2 ), (seams( i, 3 ) + 1) % 3 ]
        boundaries: Edges with only one incident triangle, as a #boundaries-by-2 matrix of
                    indices. Each row is organized as [ face_index, face_vertex_index ]
                    such that the edge is:
                        F[ boundaries( i, 0 ), boundaries( i, 1 ) ], F[ boundaries( i, 0 ), (boundaries( i, 1 ) + 1) % 3 ]
        foldovers: Edges where the two incident triangles fold over each other in UV-space,
                   as a #foldovers-by-4 matrix of indices.
                   Each row is organized as [ forward_face_index, forward_face_vertex_index,
                   backwards_face_index, backwards_face_vertex_index ]
                   such that one side of the foldover is the edge:
                        F[ foldovers( i, 0 ), foldovers( i, 1 ) ], F[ foldovers( i, 0 ), (foldovers( i, 1 ) + 1) % 3 ]
                   and the other side is the edge:
                        F[ foldovers( i, 2 ), foldovers( i, 3 ) ], F[ foldovers( i, 2 ), (foldovers( i, 3 ) + 1) % 3 ]

    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>
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;
    // 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
    const auto& Orientation = []( const Vector2S& a, const Vector2S& b, const Vector2S& c ) -> typename DerivedT::Scalar
    {
        const Vector2S row0 = a - c;
        const Vector2S row1 = b - c;
        return row0(0)*row1(1) - row1(0)*row0(1);
    };

    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();
	const 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 a canonical orientation 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 ) {
        // The canonical orientation is the one where the smaller of
        // the two vertex indices is first.
        undirected_position_edges.insert( std::make_pair( std::min( el.first.first, el.first.second ), std::max( el.first.first, el.first.second ) ) );
    }

    // 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 ) {
        // We should only see canonical edges,
        // where the first vertex index is smaller.
        assert( vp_edge.first < vp_edge.second );

        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 );

            // 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() );
}

#endif // SEAM_EDGES_H

## seam_edges_test.cpp
// Copyright (C) 2016 Yotam Gingold <yotam@yotamgingold.com>
// GitHub Gist: https://gist.github.com/yig/bea834b6e5262490611cdb4470a2f34c
// c++ -std=c++11 -I/usr/local/include/eigen3/ seam_edges_test.cpp -o seam_edges_test

/*
Boost Software License - Version 1.0 - August 17th, 2003

Permission is hereby granted, free of charge, to any person or organization
obtaining a copy of the software and accompanying documentation covered by
this license (the "Software") to use, reproduce, display, distribute,
execute, and transmit the Software, and to prepare derivative works of the
Software, and to permit third-parties to whom the Software is furnished to
do so, all subject to the following:

The copyright notices in the Software and this entire statement, including
the above license grant, this restriction and the following disclaimer,
must be included in all copies of the Software, in whole or in part, and
all derivative works of the Software, unless such copies or derivative
works are solely in the form of machine-executable object code generated by
a source language processor.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
DEALINGS IN THE SOFTWARE.
*/

#include <iostream>
#include <Eigen/Core>
#include "seam_edges.h"

using namespace std;
using namespace Eigen;

template< typename DerivedV, typename DerivedTC, typename DerivedF >
void test_mesh( const DerivedV& V, const DerivedTC& TC, const DerivedF& F, const DerivedF& FTC ) {
	MatrixXi seams, boundaries, foldovers;
	seam_edges( V, TC, F, FTC, seams, boundaries, foldovers );

	cout << "seams:\n";
	cout << seams << "\n\n";

	cout << "boundaries:\n";
	cout << boundaries << "\n\n";

	cout << "foldovers:\n";
	cout << foldovers << "\n\n";
}

void test_one_triangle() {
	MatrixXd V( 3, 3 );
	V.row(0) << 0, 0, 0;
	V.row(1) << 1, 0, 0;
	V.row(2) << 0, 1, 0;

	MatrixXd TC( 3, 2 );
	TC = V.block< 3, 2 >( 0, 0 );

	MatrixXi F( 1, 3 );
	F << 0, 1, 2;

	// Should print every edge of triangle 0 as a boundary:
	/*
    0 0
    0 1
    0 2
	*/
	test_mesh( V, TC, F, F );
}

void test_cube() {
	MatrixXd V( 8, 3 );
	V.row(0) << 0.0, 0.0, 0.0;
	V.row(1) << 0.0, 0.0, 1.0;
	V.row(2) << 0.0, 1.0, 0.0;
	V.row(3) << 0.0, 1.0, 1.0;
	V.row(4) << 1.0, 0.0, 0.0;
	V.row(5) << 1.0, 0.0, 1.0;
	V.row(6) << 1.0, 1.0, 0.0;
	V.row(7) << 1.0, 1.0, 1.0;

	MatrixXd TC( 14, 2 );
    TC.row(0)  << 0.25 , 0.75;
    TC.row(1)  << 0.5  , 0.75;
    TC.row(2)  << 0    , 0.5;
    TC.row(3)  << 0.25 , 0.5;
    TC.row(4)  << 0.5  , 0.5;
    TC.row(5)  << 0.75 , 0.5;
    TC.row(6)  << 1.0  , 0.5;
    TC.row(7)  << 0.0  , 0.25;
    TC.row(8)  << 0.25 , 0.25;
    TC.row(9)  << 0.5  , 0.25;
    TC.row(10) << 0.75 , 0.25;
    TC.row(11) << 1.0  , 0.25;
    TC.row(12) << 0.25 , 0.0;
    TC.row(13) << 0.5  , 0.0;

	MatrixXi F( 12, 3 );
    F.row(0)    << 0, 6, 4;
    F.row(1)    << 0, 2, 6;
    F.row(2)    << 0, 3, 2;
    F.row(3)    << 0, 1, 3;
    F.row(4)    << 2, 7, 6;
    F.row(5)    << 2, 3, 7;
    F.row(6)    << 4, 6, 7;
    F.row(7)    << 4, 7, 5;
    F.row(8)    << 0, 4, 5;
    F.row(9)    << 0, 5, 1;
    F.row(10)   << 1, 5, 7;
    F.row(11)   << 1, 7, 3;

	MatrixXi FTC( 12, 3 );
    FTC.row(0)  << 9  , 12 , 13;
    FTC.row(1)  << 9  , 8  , 12;
    FTC.row(2)  << 9  , 3  , 8;
    FTC.row(3)  << 9  , 4  , 3;
    FTC.row(4)  << 8  , 2  , 7;
    FTC.row(5)  << 8  , 3  , 2;
    FTC.row(6)  << 10 , 11 , 6;
    FTC.row(7)  << 10 , 6  , 5;
    FTC.row(8)  << 9  , 10 , 5;
    FTC.row(9)  << 9  , 5  , 4;
    FTC.row(10) << 4  , 1  , 0;
    FTC.row(11) << 4  , 0  , 3;

	// Should find only seams:
    /*
    1  1  4  2
    6  1  4  1
    8  0  0  2
    10 1  7  1
    10 0  9  1
    6  0  0  1
    5  1 11  1
    */
	test_mesh( V, TC, F, FTC );
}

int main(int argc, char *argv[]) {
    cout << "=> test_one_triangle()\n";
	test_one_triangle();

    cout << endl << endl;

    cout << "=> test_cube()\n";
    test_cube();
}
	// Copyright (C) 2016 Yotam Gingold <yotam@yotamgingold.com>
	// GitHub Gist: https://gist.github.com/yig/bea834b6e5262490611cdb4470a2f34c

	/*
	Boost Software License - Version 1.0 - August 17th, 2003

	Permission is hereby granted, free of charge, to any person or organization
	obtaining a copy of the software and accompanying documentation covered by
	this license (the "Software") to use, reproduce, display, distribute,
	execute, and transmit the Software, and to prepare derivative works of the
	Software, and to permit third-parties to whom the Software is furnished to
	do so, all subject to the following:

	The copyright notices in the Software and this entire statement, including
	the above license grant, this restriction and the following disclaimer,
	must be included in all copies of the Software, in whole or in part, and
	all derivative works of the Software, unless such copies or derivative
	works are solely in the form of machine-executable object code generated by
	a source language processor.

	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
	FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
	SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
	FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
	ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
	DEALINGS IN THE SOFTWARE.
	*/

	#ifndef SEAM_EDGES_H
	#define SEAM_EDGES_H

	#include <Eigen/Core>
	#include <unordered_map>
	#include <unordered_set>
	#include <cassert>

	/*
	Returns all UV-space boundaries of a mesh.
	Inputs:
	V: The positions of the input mesh.
	TC: The 2D texture coordinates of the input mesh (2 columns).
	F: A manifold-with-boundary triangle mesh, as vertex indices into `V` for the three vertices of each triangle.
	FTC: Indices into `TC` for the three vertices of each triangle.
	Outputs:
	seams: Edges where the forwards and backwards directions have different texture
	coordinates, as a #seams-by-4 matrix of indices.
	Each row is organized as [ forward_face_index, forward_face_vertex_index,
	backwards_face_index, backwards_face_vertex_index ]
	such that one side of the seam is the edge:
	F[ seams( i, 0 ), seams( i, 1 ) ], F[ seams( i, 0 ), (seams( i, 1 ) + 1) % 3 ]
	and the other side is the edge:
	F[ seams( i, 2 ), seams( i, 3 ) ], F[ seams( i, 2 ), (seams( i, 3 ) + 1) % 3 ]
	boundaries: Edges with only one incident triangle, as a #boundaries-by-2 matrix of
	indices. Each row is organized as [ face_index, face_vertex_index ]
	such that the edge is:
	F[ boundaries( i, 0 ), boundaries( i, 1 ) ], F[ boundaries( i, 0 ), (boundaries( i, 1 ) + 1) % 3 ]
	foldovers: Edges where the two incident triangles fold over each other in UV-space,
	as a #foldovers-by-4 matrix of indices.
	Each row is organized as [ forward_face_index, forward_face_vertex_index,
	backwards_face_index, backwards_face_vertex_index ]
	such that one side of the foldover is the edge:
	F[ foldovers( i, 0 ), foldovers( i, 1 ) ], F[ foldovers( i, 0 ), (foldovers( i, 1 ) + 1) % 3 ]
	and the other side is the edge:
	F[ foldovers( i, 2 ), foldovers( i, 3 ) ], F[ foldovers( i, 2 ), (foldovers( i, 3 ) + 1) % 3 ]

	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>
	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;
	// 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
	const auto& Orientation = []( const Vector2S& a, const Vector2S& b, const Vector2S& c ) -> typename DerivedT::Scalar
	{
	const Vector2S row0 = a - c;
	const Vector2S row1 = b - c;
	return row0(0)row1(1) - row1(0)row0(1);
	};

	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();
	const 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 a canonical orientation 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 ) {
	// The canonical orientation is the one where the smaller of
	// the two vertex indices is first.
	undirected_position_edges.insert( std::make_pair( std::min( el.first.first, el.first.second ), std::max( el.first.first, el.first.second ) ) );
	}

	// 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 ) {
	// We should only see canonical edges,
	// where the first vertex index is smaller.
	assert( vp_edge.first < vp_edge.second );

	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 );

	// 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() );
	}

	#endif // SEAM_EDGES_H
	// Copyright (C) 2016 Yotam Gingold <yotam@yotamgingold.com>
	// GitHub Gist: https://gist.github.com/yig/bea834b6e5262490611cdb4470a2f34c
	// c++ -std=c++11 -I/usr/local/include/eigen3/ seam_edges_test.cpp -o seam_edges_test

	/*
	Boost Software License - Version 1.0 - August 17th, 2003

	Permission is hereby granted, free of charge, to any person or organization
	obtaining a copy of the software and accompanying documentation covered by
	this license (the "Software") to use, reproduce, display, distribute,
	execute, and transmit the Software, and to prepare derivative works of the
	Software, and to permit third-parties to whom the Software is furnished to
	do so, all subject to the following:

	The copyright notices in the Software and this entire statement, including
	the above license grant, this restriction and the following disclaimer,
	must be included in all copies of the Software, in whole or in part, and
	all derivative works of the Software, unless such copies or derivative
	works are solely in the form of machine-executable object code generated by
	a source language processor.

	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
	FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
	SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
	FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
	ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
	DEALINGS IN THE SOFTWARE.
	*/

	#include <iostream>
	#include <Eigen/Core>
	#include "seam_edges.h"

	using namespace std;
	using namespace Eigen;

	template< typename DerivedV, typename DerivedTC, typename DerivedF >
	void test_mesh( const DerivedV& V, const DerivedTC& TC, const DerivedF& F, const DerivedF& FTC ) {
	MatrixXi seams, boundaries, foldovers;
	seam_edges( V, TC, F, FTC, seams, boundaries, foldovers );

	cout << "seams:\n";
	cout << seams << "\n\n";

	cout << "boundaries:\n";
	cout << boundaries << "\n\n";

	cout << "foldovers:\n";
	cout << foldovers << "\n\n";
	}

	void test_one_triangle() {
	MatrixXd V( 3, 3 );
	V.row(0) << 0, 0, 0;
	V.row(1) << 1, 0, 0;
	V.row(2) << 0, 1, 0;

	MatrixXd TC( 3, 2 );
	TC = V.block< 3, 2 >( 0, 0 );

	MatrixXi F( 1, 3 );
	F << 0, 1, 2;

	// Should print every edge of triangle 0 as a boundary:
	/*
	0 0
	0 1
	0 2
	*/
	test_mesh( V, TC, F, F );
	}

	void test_cube() {
	MatrixXd V( 8, 3 );
	V.row(0) << 0.0, 0.0, 0.0;
	V.row(1) << 0.0, 0.0, 1.0;
	V.row(2) << 0.0, 1.0, 0.0;
	V.row(3) << 0.0, 1.0, 1.0;
	V.row(4) << 1.0, 0.0, 0.0;
	V.row(5) << 1.0, 0.0, 1.0;
	V.row(6) << 1.0, 1.0, 0.0;
	V.row(7) << 1.0, 1.0, 1.0;

	MatrixXd TC( 14, 2 );
	TC.row(0) << 0.25 , 0.75;
	TC.row(1) << 0.5 , 0.75;
	TC.row(2) << 0 , 0.5;
	TC.row(3) << 0.25 , 0.5;
	TC.row(4) << 0.5 , 0.5;
	TC.row(5) << 0.75 , 0.5;
	TC.row(6) << 1.0 , 0.5;
	TC.row(7) << 0.0 , 0.25;
	TC.row(8) << 0.25 , 0.25;
	TC.row(9) << 0.5 , 0.25;
	TC.row(10) << 0.75 , 0.25;
	TC.row(11) << 1.0 , 0.25;
	TC.row(12) << 0.25 , 0.0;
	TC.row(13) << 0.5 , 0.0;

	MatrixXi F( 12, 3 );
	F.row(0) << 0, 6, 4;
	F.row(1) << 0, 2, 6;
	F.row(2) << 0, 3, 2;
	F.row(3) << 0, 1, 3;
	F.row(4) << 2, 7, 6;
	F.row(5) << 2, 3, 7;
	F.row(6) << 4, 6, 7;
	F.row(7) << 4, 7, 5;
	F.row(8) << 0, 4, 5;
	F.row(9) << 0, 5, 1;
	F.row(10) << 1, 5, 7;
	F.row(11) << 1, 7, 3;

	MatrixXi FTC( 12, 3 );
	FTC.row(0) << 9 , 12 , 13;
	FTC.row(1) << 9 , 8 , 12;
	FTC.row(2) << 9 , 3 , 8;
	FTC.row(3) << 9 , 4 , 3;
	FTC.row(4) << 8 , 2 , 7;
	FTC.row(5) << 8 , 3 , 2;
	FTC.row(6) << 10 , 11 , 6;
	FTC.row(7) << 10 , 6 , 5;
	FTC.row(8) << 9 , 10 , 5;
	FTC.row(9) << 9 , 5 , 4;
	FTC.row(10) << 4 , 1 , 0;
	FTC.row(11) << 4 , 0 , 3;

	// Should find only seams:
	/*
	1 1 4 2
	6 1 4 1
	8 0 0 2
	10 1 7 1
	10 0 9 1
	6 0 0 1
	5 1 11 1
	*/
	test_mesh( V, TC, F, FTC );
	}

	int main(int argc, char *argv[]) {
	cout << "=> test_one_triangle()\n";
	test_one_triangle();

	cout << endl << endl;

	cout << "=> test_cube()\n";
	test_cube();
	}