C++ PlainObjectBase::maxCoeff方法代码示例

IGL_INLINE void igl::adjacency_matrix(
  const Eigen::PlainObjectBase<DerivedF> & F, 
  Eigen::SparseMatrix<T>& A)
{
  using namespace std;
  using namespace Eigen;
  typedef typename DerivedF::Scalar Index;

  typedef Triplet<T> IJV;
  vector<IJV > ijv;
  ijv.reserve(F.size()*2);
  // Loop over faces
  for(int i = 0;i<F.rows();i++)
  {
    // Loop over this face
    for(int j = 0;j<F.cols();j++)
    {
      // Get indices of edge: s --> d
      Index s = F(i,j);
      Index d = F(i,(j+1)%F.cols());
      ijv.push_back(IJV(s,d,1));
      ijv.push_back(IJV(d,s,1));
    }
  }

  const Index n = F.maxCoeff()+1;
  A.resize(n,n);
  switch(F.cols())
  {
    case 3:
      A.reserve(6*(F.maxCoeff()+1));
      break;
    case 4:
      A.reserve(26*(F.maxCoeff()+1));
      break;
  }
  A.setFromTriplets(ijv.begin(),ijv.end());

  // Force all non-zeros to be one

  // Iterate over outside
  for(int k=0; k<A.outerSize(); ++k)
  {
    // Iterate over inside
    for(typename Eigen::SparseMatrix<T>::InnerIterator it (A,k); it; ++it)
    {
      assert(it.value() != 0);
      A.coeffRef(it.row(),it.col()) = 1;
    }
  }
}

IGL_INLINE std::vector<bool> igl::is_irregular_vertex(const Eigen::PlainObjectBase<DerivedV> &V, const Eigen::PlainObjectBase<DerivedF> &F)
{
  Eigen::VectorXi count = Eigen::VectorXi::Zero(F.maxCoeff());

  for(unsigned i=0; i<F.rows();++i)
  {
    for(unsigned j=0; j<F.cols();++j)
    {
      if (F(i,j) < F(i,(j+1)%F.cols())) // avoid duplicate edges
      {
        count(F(i,j  )) += 1;
        count(F(i,(j+1)%F.cols())) += 1;
      }
    }
  }

  std::vector<bool> border = is_border_vertex(V,F);

  std::vector<bool> res(count.size());

  for (unsigned i=0; i<res.size(); ++i)
    res[i] = !border[i] && count[i] != (F.cols() == 3 ? 6 : 4 );

  return res;
}

IGL_INLINE void igl::jet(
  const Eigen::PlainObjectBase<DerivedZ> & Z,
  const bool normalize,
  Eigen::PlainObjectBase<DerivedC> & C)
{
  const double min_z = (normalize?Z.minCoeff():0);
  const double max_z = (normalize?Z.maxCoeff():-1);
  return jet(Z,min_z,max_z,C);
}

IGL_INLINE void igl::components(
  const Eigen::SparseMatrix<AScalar> & A,
  Eigen::PlainObjectBase<DerivedC> & C,
  Eigen::PlainObjectBase<Derivedcounts> & counts)
{
  using namespace Eigen;
  using namespace std;
  assert(A.rows() == A.cols() && "A should be square.");
  const size_t n = A.rows();
  Array<bool,Dynamic,1> seen = Array<bool,Dynamic,1>::Zero(n,1);
  C.resize(n,1);
  typename DerivedC::Scalar id = 0;
  vector<typename Derivedcounts::Scalar> vcounts;
  // breadth first search
  for(int k=0; k<A.outerSize(); ++k)
  {
    if(seen(k))
    {
      continue;
    }
    queue<int> Q;
    Q.push(k);
    vcounts.push_back(0);
    while(!Q.empty())
    {
      const int f = Q.front();
      Q.pop();
      if(seen(f))
      {
        continue;
      }
      seen(f) = true;
      C(f,0) = id;
      vcounts[id]++;
      // Iterate over inside
      for(typename SparseMatrix<AScalar>::InnerIterator it (A,f); it; ++it)
      {
        const int g = it.index();
        if(!seen(g) && it.value())
        {
          Q.push(g);
        }
      }
    }
    id++;
  }
  assert((size_t) id == vcounts.size());
  const size_t ncc = vcounts.size();
  assert((size_t)C.maxCoeff()+1 == ncc);
  counts.resize(ncc,1);
  for(size_t i = 0;i<ncc;i++)
  {
    counts(i) = vcounts[i];
  }
}

IGL_INLINE void igl::directed_edge_parents(
  const Eigen::PlainObjectBase<DerivedE> & E,
  Eigen::PlainObjectBase<DerivedP> & P)
{
  using namespace Eigen;
  using namespace std;
  VectorXi I = VectorXi::Constant(E.maxCoeff()+1,1,-1);
  //I(E.col(1)) = 0:E.rows()-1
  slice_into(colon<int>(0,E.rows()-1),E.col(1).eval(),I);
  VectorXi roots,_;
  setdiff(E.col(0).eval(),E.col(1).eval(),roots,_);
  for_each(roots.data(),roots.data()+roots.size(),[&](int r){I(r)=-1;});
  slice(I,E.col(0).eval(),P);
}

IGL_INLINE void igl::slice(
  const Eigen::PlainObjectBase<DerivedX> & X,
  const Eigen::PlainObjectBase<DerivedR> & R,
  const Eigen::PlainObjectBase<DerivedC> & C,
  Eigen::PlainObjectBase<DerivedY> & Y)
{
#ifndef NDEBUG
  int xm = X.rows();
  int xn = X.cols();
#endif
  int ym = R.size();
  int yn = C.size();

  // special case when R or C is empty
  if(ym == 0 || yn == 0)
  {
    Y.resize(ym,yn);
    return;
  }

  assert(R.minCoeff() >= 0);
  assert(R.maxCoeff() < xm);
  assert(C.minCoeff() >= 0);
  assert(C.maxCoeff() < xn);

  // Resize output
  Y.resize(ym,yn);
  // loop over output rows, then columns
  for(int i = 0;i<ym;i++)
  {
    for(int j = 0;j<yn;j++)
    {
      Y(i,j) = X(R(i),C(j));
    }
  }
}

IGL_INLINE void igl::cut_mesh(
                                                                  const Eigen::PlainObjectBase<DerivedV> &V,
                                                                  const Eigen::PlainObjectBase<DerivedF> &F,
                                                                  const std::vector<std::vector<VFType> >& VF,
                                                                  const std::vector<std::vector<VFType> >& VFi,
                                                                  const Eigen::PlainObjectBase<DerivedTT>& TT,
                                                                  const Eigen::PlainObjectBase<DerivedTT>& TTi,
                                                                  const std::vector<bool> &V_border,
                                                                  const Eigen::PlainObjectBase<DerivedC> &cuts,
                                                                  Eigen::PlainObjectBase<DerivedV> &Vcut,
                                                                  Eigen::PlainObjectBase<DerivedF> &Fcut)
{
  //finding the cuts is done, now we need to actually generate a cut mesh
  igl::MeshCutterMini<DerivedV, DerivedF, VFType, DerivedTT, DerivedC> mc(V, F, TT, TTi, VF, VFi, V_border, cuts);
  mc.InitMappingSeam();

  Fcut = mc.HandleS_Index;
  //we have the faces, we need the vertices;
  int newNumV = Fcut.maxCoeff()+1;
  Vcut.setZero(newNumV,3);
  for (int vi=0; vi<V.rows(); ++vi)
    for (int i=0; i<mc.HandleV_Integer[vi].size();++i)
      Vcut.row(mc.HandleV_Integer[vi][i]) = V.row(vi);

  //ugly hack to fix some problematic cases (border vertex that is also on the boundary of the hole
  for (int fi =0; fi<Fcut.rows(); ++fi)
    for (int k=0; k<3; ++k)
      if (Fcut(fi,k)==-1)
      {
        //we need to add a vertex
        Fcut(fi,k) = newNumV;
        newNumV ++;
        Vcut.conservativeResize(newNumV, Eigen::NoChange);
        Vcut.row(newNumV-1) = V.row(F(fi,k));
      }


}

IGL_INLINE void igl::gaussian_curvature(
  const Eigen::PlainObjectBase<DerivedV>& V,
  const Eigen::PlainObjectBase<DerivedF>& F,
  Eigen::PlainObjectBase<DerivedK> & K)
{
  using namespace Eigen;
  using namespace std;
  // internal corner angles
  Matrix<
    typename DerivedV::Scalar,
    DerivedF::RowsAtCompileTime,
    DerivedF::ColsAtCompileTime> A;
  internal_angles(V,F,A);
  K.resize(V.rows(),1);
  K.setConstant(V.rows(),1,2.*PI);
  assert(A.rows() == F.rows());
  assert(A.cols() == F.cols());
  assert(K.rows() == V.rows());
  assert(F.maxCoeff() < V.rows());
  assert(K.cols() == 1);
  const int Frows = F.rows();
  //K_G(x_i) = (2π - ∑θj)
//#ifndef IGL_GAUSSIAN_CURVATURE_OMP_MIN_VALUE
//#  define IGL_GAUSSIAN_CURVATURE_OMP_MIN_VALUE 1000
//#endif
//#pragma omp parallel for if (Frows>IGL_GAUSSIAN_CURVATURE_OMP_MIN_VALUE)
  for(int f = 0;f<Frows;f++)
  {
    // throw normal at each corner
    for(int j = 0; j < 3;j++)
    {
      // Q: Does this need to be critical?
      // H: I think so, sadly. Maybe there's a way to use reduction
//#pragma omp critical
      K(F(f,j),0) -=  A(f,j);
    }
  }
}

IGL_INLINE bool igl::arap_precomputation(
    const Eigen::PlainObjectBase<DerivedV> & V,
    const Eigen::PlainObjectBase<DerivedF> & F,
    const int dim,
    const Eigen::PlainObjectBase<Derivedb> & b,
    ARAPData & data)
{
    using namespace std;
    using namespace Eigen;
    typedef typename DerivedV::Scalar Scalar;
    // number of vertices
    const int n = V.rows();
    data.n = n;
    assert((b.size() == 0 || b.maxCoeff() < n) && "b out of bounds");
    assert((b.size() == 0 || b.minCoeff() >=0) && "b out of bounds");
    // remember b
    data.b = b;
    //assert(F.cols() == 3 && "For now only triangles");
    // dimension
    //const int dim = V.cols();
    assert((dim == 3 || dim ==2) && "dim should be 2 or 3");
    data.dim = dim;
    //assert(dim == 3 && "Only 3d supported");
    // Defaults
    data.f_ext = MatrixXd::Zero(n,data.dim);

    assert(data.dim <= V.cols() && "solve dim should be <= embedding");
    bool flat = (V.cols() - data.dim)==1;

    DerivedV plane_V;
    DerivedF plane_F;
    typedef SparseMatrix<Scalar> SparseMatrixS;
    SparseMatrixS ref_map,ref_map_dim;
    if(flat)
    {
        project_isometrically_to_plane(V,F,plane_V,plane_F,ref_map);
        repdiag(ref_map,dim,ref_map_dim);
    }
    const PlainObjectBase<DerivedV>& ref_V = (flat?plane_V:V);
    const PlainObjectBase<DerivedF>& ref_F = (flat?plane_F:F);
    SparseMatrixS L;
    cotmatrix(V,F,L);

    ARAPEnergyType eff_energy = data.energy;
    if(eff_energy == ARAP_ENERGY_TYPE_DEFAULT)
    {
        switch(F.cols())
        {
        case 3:
            if(data.dim == 3)
            {
                eff_energy = ARAP_ENERGY_TYPE_SPOKES_AND_RIMS;
            } else
            {
                eff_energy = ARAP_ENERGY_TYPE_ELEMENTS;
            }
            break;
        case 4:
            eff_energy = ARAP_ENERGY_TYPE_ELEMENTS;
            break;
        default:
            assert(false);
        }
    }


    // Get covariance scatter matrix, when applied collects the covariance
    // matrices used to fit rotations to during optimization
    covariance_scatter_matrix(ref_V,ref_F,eff_energy,data.CSM);
    if(flat)
    {
        data.CSM = (data.CSM * ref_map_dim.transpose()).eval();
    }
    assert(data.CSM.cols() == V.rows()*data.dim);

    // Get group sum scatter matrix, when applied sums all entries of the same
    // group according to G
    SparseMatrix<double> G_sum;
    if(data.G.size() == 0)
    {
        if(eff_energy == ARAP_ENERGY_TYPE_ELEMENTS)
        {
            speye(F.rows(),G_sum);
        } else
        {
            speye(n,G_sum);
        }
    } else
    {
        // groups are defined per vertex, convert to per face using mode
        if(eff_energy == ARAP_ENERGY_TYPE_ELEMENTS)
        {
            Eigen::Matrix<int,Eigen::Dynamic,1> GG;
            MatrixXi GF(F.rows(),F.cols());
            for(int j = 0; j<F.cols(); j++)
            {
                Matrix<int,Eigen::Dynamic,1> GFj;
                slice(data.G,F.col(j),GFj);
                GF.col(j) = GFj;
            }
//.........这里部分代码省略.........

IGL_INLINE void igl::adjacency_list(
    const Eigen::PlainObjectBase<Index>  & F,
    std::vector<std::vector<IndexVector> >& A,
    bool sorted)
{
  A.clear(); 
  A.resize(F.maxCoeff()+1);
  
  // Loop over faces
  for(int i = 0;i<F.rows();i++)
  {
    // Loop over this face
    for(int j = 0;j<F.cols();j++)
    {
      // Get indices of edge: s --> d
      int s = F(i,j);
      int d = F(i,(j+1)%F.cols());
      A.at(s).push_back(d);
      A.at(d).push_back(s);
    }
  }
  
  // Remove duplicates
  for(int i=0; i<(int)A.size();++i)
  {
    std::sort(A[i].begin(), A[i].end());
    A[i].erase(std::unique(A[i].begin(), A[i].end()), A[i].end());
  }
  
  // If needed, sort every VV
  if (sorted)
  {
    // Loop over faces
    
    // for every vertex v store a set of ordered edges not incident to v that belongs to triangle incident on v.
    std::vector<std::vector<std::vector<int> > > SR; 
    SR.resize(A.size());
    
    for(int i = 0;i<F.rows();i++)
    {
      // Loop over this face
      for(int j = 0;j<F.cols();j++)
      {
        // Get indices of edge: s --> d
        int s = F(i,j);
        int d = F(i,(j+1)%F.cols());
        // Get index of opposing vertex v
        int v = F(i,(j+2)%F.cols());
        
        std::vector<int> e(2);
        e[0] = d;
        e[1] = v;
        SR[s].push_back(e);
      }
    }
    
    for(int v=0; v<(int)SR.size();++v)
    {
      std::vector<IndexVector>& vv = A.at(v);
      std::vector<std::vector<int> >& sr = SR[v];
      
      std::vector<std::vector<int> > pn = sr;
      
      // Compute previous/next for every element in sr
      for(int i=0;i<(int)sr.size();++i)
      {
        int a = sr[i][0];
        int b = sr[i][1];
        
        // search for previous
        int p = -1;
        for(int j=0;j<(int)sr.size();++j)
          if(sr[j][1] == a)
            p = j;
        pn[i][0] = p;
        
        // search for next
        int n = -1;
        for(int j=0;j<(int)sr.size();++j)
          if(sr[j][0] == b)
            n = j;
        pn[i][1] = n;
        
      }
      
      // assume manifoldness (look for beginning of a single chain)
      int c = 0;
      for(int j=0; j<=(int)sr.size();++j)
        if (pn[c][0] != -1)
          c = pn[c][0];
      
      if (pn[c][0] == -1) // border case
      {
        // finally produce the new vv relation
        for(int j=0; j<(int)sr.size();++j)
        {
          vv[j] = sr[c][0];
          if (pn[c][1] != -1)
            c = pn[c][1];
        }
//.........这里部分代码省略.........

void igl::bfs_orient(
  const Eigen::PlainObjectBase<DerivedF> & F,
  Eigen::PlainObjectBase<DerivedFF> & FF,
  Eigen::PlainObjectBase<DerivedC> & C)
{
  using namespace Eigen;
  using namespace igl;
  using namespace std;
  SparseMatrix<int> A;
  manifold_patches(F,C,A);

  // number of faces
  const int m = F.rows();
  // number of patches
  const int num_cc = C.maxCoeff()+1;
  VectorXi seen = VectorXi::Zero(m);

  // Edge sets
  const int ES[3][2] = {{1,2},{2,0},{0,1}};

  if(&FF != &F)
  {
    FF = F;
  }
  // loop over patches
#pragma omp parallel for
  for(int c = 0;c<num_cc;c++)
  {
    queue<int> Q;
    // find first member of patch c
    for(int f = 0;f<FF.rows();f++)
    {
      if(C(f) == c)
      {
        Q.push(f);
        break;
      }
    }
    assert(!Q.empty());
    while(!Q.empty())
    {
      const int f = Q.front();
      Q.pop();
      if(seen(f) > 0)
      {
        continue;
      }
      seen(f)++;
      // loop over neighbors of f
      for(typename SparseMatrix<int>::InnerIterator it (A,f); it; ++it)
      {
        // might be some lingering zeros, and skip self-adjacency
        if(it.value() != 0 && it.row() != f)
        {
          const int n = it.row();
          assert(n != f);
          // loop over edges of f
          for(int efi = 0;efi<3;efi++)
          {
            // efi'th edge of face f
            Vector2i ef(FF(f,ES[efi][0]),FF(f,ES[efi][1]));
            // loop over edges of n
            for(int eni = 0;eni<3;eni++)
            {
              // eni'th edge of face n
              Vector2i en(FF(n,ES[eni][0]),FF(n,ES[eni][1]));
              // Match (half-edges go same direction)
              if(ef(0) == en(0) && ef(1) == en(1))
              {
                // flip face n
                FF.row(n) = FF.row(n).reverse().eval();
              }
            }
          }
          // add neighbor to queue
          Q.push(n);
        }
      }
    }
  }

  // make sure flip is OK if &FF = &F
}

IGL_INLINE void igl::reorient_facets_raycast(
  const Eigen::PlainObjectBase<DerivedV> & V,
  const Eigen::PlainObjectBase<DerivedF> & F,
  int rays_total,
  int rays_minimum,
  bool facet_wise,
  bool use_parity,
  bool is_verbose,
  Eigen::PlainObjectBase<DerivedI> & I,
  Eigen::PlainObjectBase<DerivedC> & C)
{
  using namespace Eigen;
  using namespace std;
  assert(F.cols() == 3);
  assert(V.cols() == 3);
  
  // number of faces
  const int m = F.rows();
  
  MatrixXi FF = F;
  if (facet_wise) {
    C.resize(m);
    for (int i = 0; i < m; ++i) C(i) = i;
  
  } else {
    if (is_verbose) cout << "extracting patches... ";
    bfs_orient(F,FF,C);
  }
  if (is_verbose) cout << (C.maxCoeff() + 1)  << " components. ";
  
  // number of patches
  const int num_cc = C.maxCoeff()+1;
  
  // Init Embree
  EmbreeIntersector ei;
  ei.init(V.template cast<float>(),FF);
  
  // face normal
  MatrixXd N;
  per_face_normals(V,FF,N);
  
  // face area
  Matrix<typename DerivedV::Scalar,Dynamic,1> A;
  doublearea(V,FF,A);
  double area_total = A.sum();
  
  // determine number of rays per component according to its area
  VectorXd area_per_component;
  area_per_component.setZero(num_cc);
  for (int f = 0; f < m; ++f)
  {
    area_per_component(C(f)) += A(f);
  }
  VectorXi num_rays_per_component(num_cc);
  for (int c = 0; c < num_cc; ++c)
  {
    num_rays_per_component(c) = max<int>(static_cast<int>(rays_total * area_per_component(c) / area_total), rays_minimum);
  }
  rays_total = num_rays_per_component.sum();
  
  // generate all the rays
  if (is_verbose) cout << "generating rays... ";
  uniform_real_distribution<float> rdist;
  mt19937 prng;
  prng.seed(time(nullptr));
  vector<int     > ray_face;
  vector<Vector3f> ray_ori;
  vector<Vector3f> ray_dir;
  ray_face.reserve(rays_total);
  ray_ori .reserve(rays_total);
  ray_dir .reserve(rays_total);
  for (int c = 0; c < num_cc; ++c)
  {
    if (area_per_component[c] == 0)
    {
      continue;
    }
    vector<int> CF;     // set of faces per component
    vector<double> CF_area;
    for (int f = 0; f < m; ++f)
    {
      if (C(f)==c)
      {
        CF.push_back(f);
        CF_area.push_back(A(f));
      }
    }
    // discrete distribution for random selection of faces with probability proportional to their areas
    discrete_distribution<int> ddist(CF.size(), 0, CF.size(), [&](double i){ return CF_area[static_cast<int>(i)]; });       // simple ctor of (Iter, Iter) not provided by the stupid VC11/12
    for (int i = 0; i < num_rays_per_component[c]; ++i)
    {
      int f = CF[ddist(prng)];          // select face with probability proportional to face area
      float s = rdist(prng);            // random barycentric coordinate (reference: Generating Random Points in Triangles [Turk, Graphics Gems I 1990])
      float t = rdist(prng);
      float sqrt_t = sqrtf(t);
      float a = 1 - sqrt_t;
      float b = (1 - s) * sqrt_t;
      float c = s * sqrt_t;
      Vector3f p = a * V.row(FF(f,0)).template cast<float>().eval()       // be careful with the index!!!
                 + b * V.row(FF(f,1)).template cast<float>().eval()
//.........这里部分代码省略.........

IGL_INLINE void igl::boundary_loop(
    const Eigen::PlainObjectBase<DerivedF> & F, 
    std::vector<std::vector<Index> >& L)
{
  using namespace std;
  using namespace Eigen;
  using namespace igl;

  MatrixXd Vdummy(F.maxCoeff(),1);
  MatrixXi TT,TTi;
  vector<std::vector<int> > VF, VFi;
  triangle_triangle_adjacency(Vdummy,F,TT,TTi);
  vertex_triangle_adjacency(Vdummy,F,VF,VFi);

  vector<bool> unvisited = is_border_vertex(Vdummy,F);
  set<int> unseen;
  for (int i = 0; i < unvisited.size(); ++i)
  {
    if (unvisited[i])
      unseen.insert(unseen.end(),i);
  }

  while (!unseen.empty())
  {
    vector<Index> l;

    // Get first vertex of loop
    int start = *unseen.begin();
    unseen.erase(unseen.begin());
    unvisited[start] = false;
    l.push_back(start);

    bool done = false;
    while (!done)
    {
      // Find next vertex
      bool newBndEdge = false;
      int v = l[l.size()-1];
      int next;
      for (int i = 0; i < (int)VF[v].size() && !newBndEdge; i++)
      {
        int fid = VF[v][i];

        if (TT.row(fid).minCoeff() < 0.) // Face contains boundary edge
        {
          int vLoc;
          if (F(fid,0) == v) vLoc = 0;
          if (F(fid,1) == v) vLoc = 1;
          if (F(fid,2) == v) vLoc = 2;

          int vPrev = F(fid,(vLoc + F.cols()-1) % F.cols());
          int vNext = F(fid,(vLoc + 1) % F.cols());

          bool newBndEdge = false;
          if (unvisited[vPrev] && TT(fid,(vLoc+2) % F.cols()) < 0)
          {
            next = vPrev;
            newBndEdge = true;
          }
          else if (unvisited[vNext] && TT(fid,vLoc) < 0)
          {
            next = vNext;
            newBndEdge = true;
          }
        }
      }

      if (newBndEdge)
      {
        l.push_back(next);
        unseen.erase(next);
        unvisited[next] = false;
      }
      else
        done = true;
    }
    L.push_back(l);
  }
}

IGL_INLINE void igl::orient_outward(
  const Eigen::PlainObjectBase<DerivedV> & V,
  const Eigen::PlainObjectBase<DerivedF> & F,
  const Eigen::PlainObjectBase<DerivedC> & C,
  Eigen::PlainObjectBase<DerivedFF> & FF,
  Eigen::PlainObjectBase<DerivedI> & I)
{
  using namespace Eigen;
  using namespace std;
  assert(C.rows() == F.rows());
  assert(F.cols() == 3);
  assert(V.cols() == 3);

  // number of faces
  const int m = F.rows();
  // number of patches
  const int num_cc = C.maxCoeff()+1;
  I.resize(num_cc);
  if(&FF != &F)
  {
    FF = F;
  }
  PlainObjectBase<DerivedV> N,BC,BCmean;
  Matrix<typename DerivedV::Scalar,Dynamic,1> A;
  VectorXd totA(num_cc), dot(num_cc);
  Matrix<typename DerivedV::Scalar,3,1> Z(1,1,1);
  per_face_normals(V,F,Z.normalized(),N);
  barycenter(V,F,BC);
  doublearea(V,F,A);
  BCmean.setConstant(num_cc,3,0);
  dot.setConstant(num_cc,1,0);
  totA.setConstant(num_cc,1,0);
  // loop over faces
  for(int f = 0;f<m;f++)
  {
    BCmean.row(C(f)) += A(f)*BC.row(f);
    totA(C(f))+=A(f);
  }
  // take area weighted average
  for(int c = 0;c<num_cc;c++)
  {
    BCmean.row(c) /= (typename DerivedV::Scalar) totA(c);
  }
  // subtract bcmean
  for(int f = 0;f<m;f++)
  {
    BC.row(f) -= BCmean.row(C(f));
    dot(C(f)) += A(f)*N.row(f).dot(BC.row(f));
  }
  // take area weighted average
  for(int c = 0;c<num_cc;c++)
  {
    dot(c) /= (typename DerivedV::Scalar) totA(c);
    if(dot(c) < 0)
    {
      I(c) = true;
    }else
    {
      I(c) = false;
    }
  }
  // flip according to I
  for(int f = 0;f<m;f++)
  {
    if(I(C(f)))
    {
      FF.row(f) = FF.row(f).reverse().eval();
    }
  }
}