libigl/include/igl/isolines_intrinsic.cpp

// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2023 Alec Jacobson <[email protected]>
//
// 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 "isolines_intrinsic.h"
#include "edge_crossings.h"
#include "cat.h"
#include "unique_edge_map.h"
#ifndef NDEBUG
#  include "is_edge_manifold.h"
#  include "is_vertex_manifold.h"
#endif
#include <unordered_map>
#include <vector>

template <
  typename DerivedF,
  typename DerivedS,
  typename Derivedvals,
  typename DerivediB,
  typename DerivediFI,
  typename DerivediE,
  typename DerivedI>
void igl::isolines_intrinsic(
  const Eigen::MatrixBase<DerivedF> & F,
  const Eigen::MatrixBase<DerivedS> & S,
  const Eigen::MatrixBase<Derivedvals> & vals,
  Eigen::PlainObjectBase<DerivediB> & iB,
  Eigen::PlainObjectBase<DerivediFI> & iFI,
  Eigen::PlainObjectBase<DerivediE> & iE,
  Eigen::PlainObjectBase<DerivedI> & I)
{
  using Index = typename DerivedF::Scalar;
  Eigen::Matrix<Index,Eigen::Dynamic,Eigen::Dynamic> uE;
  Eigen::Vector<Index,Eigen::Dynamic> EMAP,uEC,uEE;

  {
    Eigen::Matrix<Index,Eigen::Dynamic,Eigen::Dynamic> E;
    igl::unique_edge_map(F,E,uE,EMAP,uEC,uEE);
  }

  {
    std::vector<DerivediB> viB(vals.size());
    std::vector<DerivediFI> viFI(vals.size());
    std::vector<DerivediE> viE(vals.size());
    std::vector<DerivedI> vI(vals.size());
    int num_vertices = 0;
    for(int j = 0;j<vals.size();j++)
    {
      isolines_intrinsic(F,S,uE,EMAP,uEC,uEE,vals(j),viB[j],viFI[j],viE[j]);
      viE[j].array() += num_vertices;
      num_vertices += viB[j].rows();
      vI[j] = DerivedI::Constant(viE[j].rows(),j);
    }
    igl::cat(1,viB,iB);
    igl::cat(1,viFI,iFI);
    igl::cat(1,viE,iE);
    igl::cat(1,vI,I);
  }
}

template <
  typename DerivedF,
  typename DerivedS,
  typename DeriveduE,
  typename DerivedEMAP,
  typename DeriveduEC,
  typename DeriveduEE,
  typename DerivediB,
  typename DerivediFI,
  typename DerivediE>
void igl::isolines_intrinsic(
  const Eigen::MatrixBase<DerivedF> & F,
  const Eigen::MatrixBase<DerivedS> & S,
  const Eigen::MatrixBase<DeriveduE> & uE,
  const Eigen::MatrixBase<DerivedEMAP> & EMAP,
  const Eigen::MatrixBase<DeriveduEC> & uEC,
  const Eigen::MatrixBase<DeriveduEE> & uEE,
  const typename DerivedS::Scalar val,
  Eigen::PlainObjectBase<DerivediB> & iB,
  Eigen::PlainObjectBase<DerivediFI> & iFI,
  Eigen::PlainObjectBase<DerivediE> & iE)
{
  using Scalar = typename DerivedS::Scalar;

  std::unordered_map<int,int> uE2I;
  Eigen::Matrix<Scalar,Eigen::Dynamic,1> T;
  igl::edge_crossings(uE,S,val,uE2I,T);

  iB.resize(uE2I.size(),F.cols());
  iFI.resize(uE2I.size());
  Eigen::VectorXi U(uE2I.size());
  for(auto & pair : uE2I)
  {
    const int u = pair.first;
    const int w = pair.second;
    // first face incident on uE(u,:)
    const int e = uEE(uEC(u));
    const int f = e % F.rows();
    const int k = e / F.rows();
    const bool flip = uE(u,0) != F(f,(k+1)%3);
    const double t = T(w);
    iB(w,k) = 0;
    iB(w,(k+1)%3) = flip?  t:1-t;
    iB(w,(k+2)%3) = flip?1-t:t;
    iFI(w) = f;
    U(w) = u;
  }
  
  
  // Vertex crossings
  std::unordered_map<int,int> V2I;
  {
    const auto add_vertex_crossing = [&iB,&iFI](const int k, const int i, const int j)
    {
      if(k >= iB.rows())
      {
        iB.conservativeResize(2*iB.rows()+1,Eigen::NoChange);
        iFI.conservativeResize(2*iFI.rows()+1,Eigen::NoChange);
      }
      iFI(k) = i;
      iB.row(k) << 0,0,0;
      iB(k,j) = 1;
    };
    int k = iB.rows();
    for(int i = 0;i<F.rows();i++)
    {
      for(int j = 0;j<3;j++)
      {
        const int v = F(i,j);
        if(S(v) == val)
        {
          if(V2I.find(v) == V2I.end())
          {
            V2I[v] = k;
            add_vertex_crossing(k++,i,j);
          }
        }
      }
    }
    iB.conservativeResize(k,Eigen::NoChange);
    iFI.conservativeResize(k,Eigen::NoChange);
  }

  iE.resize(uE2I.size(),2);
  const auto set_row = [&iE](const int k, const int i, const int j)
  {
    if(k >= iE.rows())
    {
      iE.conservativeResize(2*iE.rows()+1,Eigen::NoChange);
    }
    iE.row(k) << i,j;
  };
  {
    int r = 0;
    for(int f = 0;f < F.rows();f++)
    {
      // find first crossing edge
      int i;
      for(i = 0;i<3;i++)
      {
        if(uE2I.find(EMAP(f+F.rows()*i)) != uE2I.end())
        {
          break;
        }
      }
      int j;
      for(j = i+1;j<3;j++)
      {
        if(uE2I.find(EMAP(f+F.rows()*j)) != uE2I.end())
        {
          break;
        }
      }
      if(j<3)
      {
        // Connect two edge crossings.

        // other vertex
        const int k = 3-i-j;
        const int wi = uE2I[EMAP(f+F.rows()*i)];
        const int wj = uE2I[EMAP(f+F.rows()*j)];
        // flip orientation based on triangle gradient
        Scalar SFfk = S(F(f,k));
        bool flip = SFfk < val;
        flip = k%2? !flip:flip;
        if(flip)
        {
          set_row(r++,wi,wj);
        }else
        {
          set_row(r++,wj,wi);
        }
      }else if(i<3)
      {
        // The only valid vertex crossing is the opposite vertex
        const int v = F(f,i);
        // Is it a crossing?
        assert(V2I.find(v) != V2I.end());
        //if(V2I.find(v) != V2I.end())
        {
          const int wv = V2I[v];
          const int wi = uE2I[EMAP(f+F.rows()*i)];
          const bool flip = S(F(f,(i+1)%3)) > val;
          if(flip)
          {
            set_row(r++,wi,wv);
          }else
          {
            set_row(r++,wv,wi);
          }
        }
      }else
      {
        // Could have 2 vertex crossings. We're only interested if there're exactly two and if the other vertex is "above".
        int i = 0;
        for(i = 0;i<3;i++)
        {
          if(S(F(f,i)) == val)
          {
            break;
          }
        }
        int j;
        for(j = i+1;j<3;j++)
        {
          if(S(F(f,j)) == val)
          {
            break;
          }
        }
        if(j<3)
        {
          // check if the third is a crossing.
          const int k = 3-i-j;
          // Triangle is constant on the val. Skip.
          if(S(F(f,k)) == val){ continue; }
          // Is this a boundary edge?
          const int u = EMAP(f+F.rows()*k);
          const int count = uEC(u+1)-uEC(u);
          if( count == 1 || S(F(f,k)) > val)
          {
            const int wi = V2I[F(f,i)];
            const int wj = V2I[F(f,j)];
            bool flip = S(F(f,k)) < val;
            flip = k%2 ? !flip:flip;
            if(flip)
            {
              set_row(r++,wj,wi);
            }else
            {
              set_row(r++,wi,wj);
            }
          }
        }
      }
    }
    iE.conservativeResize(r,Eigen::NoChange);
  }

#ifdef IGL_ISOLINES_INTRINSIC_DEBUG
  if(igl::is_vertex_manifold(F) && igl::is_edge_manifold(F))
  {
    // Check that every vertex has one in one out
    Eigen::VectorXi in_count = Eigen::VectorXi::Zero(iB.rows());
    Eigen::VectorXi out_count = Eigen::VectorXi::Zero(iB.rows());
    for(int e = 0;e<iE.rows();e++)
    {
      const int i = iE(e,0);
      out_count(i)++;
      const int j = iE(e,1);
      in_count(j)++;
    }
    for(int i = 0;i<iB.rows();i++)
    {
      assert(in_count(i) <= 1);
      assert(out_count(i) <= 1);
    }
  }
#endif
}

#ifdef IGL_STATIC_LIBRARY
// Explicit template instantiation
// generated by autoexplicit.sh
template void igl::isolines_intrinsic<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
#endif