vdpmanalyzer.cc

penMesh is a generic and efficient data structure for representing and manipulating polygonal meshes

    lp = mesh_.point(vhierarchy_.vertex_handle(*n_it));
    max_distance = std::max(max_distance, (p - lp).length());
  }

  vhierarchy_.node(node_handle).set_radius(max_distance);
}


// ----------------------------------------------------------------------------

void
compute_cone_of_normals(VHierarchyNodeHandle           node_handle,
			VHierarchyNodeHandleContainer &leaf_nodes)
{
  float           max_angle, angle;
  Vec3f           n, ln;
  VertexHandle    vh = vhierarchy_.node(node_handle).vertex_handle();
  VHierarchyNodeHandleContainer::iterator  n_it, n_end(leaf_nodes.end());

  for (Mesh::VertexFaceIter vf_it=mesh_.vf_iter(vh); vf_it; ++vf_it)
  {
    n = mesh_.calc_face_normal(vf_it);
    mesh_.set_normal(vf_it, n);
  }

  n         = mesh_.calc_vertex_normal(vh);
  max_angle = 0.0f;

  n_it = leaf_nodes.begin();
  while( n_it != n_end )
  {
    ln        = vhierarchy_.node(*n_it).normal();
    angle     = acosf( dot(n,ln) );
    max_angle = std::max(max_angle, angle );

    ++n_it;
  }

  max_angle = std::min(max_angle, float(M_PI_2));
  mesh_.set_normal(vh, n);
  vhierarchy_.node(node_handle).set_normal(n);
  vhierarchy_.node(node_handle).set_semi_angle(max_angle);
}


// ----------------------------------------------------------------------------

void
compute_screen_space_error(VHierarchyNodeHandle node_handle, VHierarchyNodeHandleContainer &leaf_nodes)
{
  std::vector<Vec3f>    residuals;
  Mesh::VertexFaceIter  vf_it;
  Mesh::HalfedgeHandle  heh;
  Mesh::VertexHandle    vh;
  Vec3f                 residual, res;
  Vec3f                 lp, tri[3];
  float                 min_distance;
  float                 s, t;
  VHierarchyNodeHandleContainer::iterator  n_it, n_end(leaf_nodes.end());

  for ( n_it = leaf_nodes.begin(); n_it != n_end; ++n_it )
  {
    lp = mesh_.point(vhierarchy_.node(*n_it).vertex_handle());

    // compute residual of a leaf-vertex from the current mesh_
    vh = vhierarchy_.node(node_handle).vertex_handle();
    residual = lp - mesh_.point(vh);
    min_distance = residual.length();

    for (vf_it=mesh_.vf_iter(vh); vf_it; ++vf_it)
    {
      heh    = mesh_.halfedge_handle(vf_it.handle());
      tri[0] = mesh_.point(mesh_.to_vertex_handle(heh));
      heh    = mesh_.next_halfedge_handle(heh);
      tri[1] = mesh_.point(mesh_.to_vertex_handle(heh));
      heh    = mesh_.next_halfedge_handle(heh);
      tri[2] = mesh_.point(mesh_.to_vertex_handle(heh));

      res = point2triangle_residual(lp, tri, s, t);

      if (res.length() < min_distance)
			{
				residual = res;
				min_distance = res.length();
			}
    }

    residuals.push_back(residual);
  }

  compute_mue_sigma(node_handle, residuals);
}


// ----------------------------------------------------------------------------

void
compute_mue_sigma(VHierarchyNodeHandle node_handle,
		  ResidualContainer &residuals)
{
  Vec3f   vn;
  float             max_inner, max_cross;
  ResidualContainer::iterator  r_it, r_end(residuals.end());

  max_inner = max_cross = 0.0f;
  vn = mesh_.normal(vhierarchy_.node(node_handle).vertex_handle());
  for (r_it = residuals.begin(); r_it != r_end; ++r_it)
  {
    float inner = fabsf(dot(*r_it, vn));
    float cross = OpenMesh::cross(*r_it, vn).length();

    max_inner = std::max(max_inner, inner);
    max_cross = std::max(max_cross, cross);
  }

  if (max_cross < 1.0e-7)
  {
    vhierarchy_.node(node_handle).set_mue(max_cross);
    vhierarchy_.node(node_handle).set_sigma(max_inner);
  }
  else {
    float  ratio = std::max(1.0f, max_inner/max_cross);
    float  whole_degree = acosf(1.0f/ratio);
    float  mue, max_mue;
    float  degree;
    float  res_length;
    Vec3f  res;

    max_mue = 0.0f;
    for (r_it = residuals.begin(); r_it != r_end; ++r_it)
    {
      res = *r_it;
      res_length = res.length();

      // TODO: take care when res.length() is too small
      degree = acosf(dot(vn,res) / res_length);

      if (degree < 0.0f)    degree = -degree;
      if (degree > float(M_PI_2))  degree = float(M_PI) - degree;

      if (degree < whole_degree)
        mue = cosf(whole_degree - degree) * res_length;
      else
        mue = res_length;

      max_mue = std::max(max_mue, mue);
    }

    vhierarchy_.node(node_handle).set_mue(max_mue);
    vhierarchy_.node(node_handle).set_sigma(ratio*max_mue);
  }
}

// ----------------------------------------------------------------------------


Vec3f
point2triangle_residual(const Vec3f &p, const Vec3f tri[3], float &s, float &t)
{ 
  OpenMesh::Vec3f B = tri[0];             // Tri.Origin();
  OpenMesh::Vec3f E0 = tri[1] - tri[0];   // rkTri.Edge0()
  OpenMesh::Vec3f E1 = tri[2] - tri[0];   // rkTri.Edge1()
  OpenMesh::Vec3f D = tri[0] - p;         // kDiff
  float	a = dot(E0, E0);                  // fA00
  float	b = dot(E0, E1);                  // fA01
  float	c = dot(E1, E1);                  // fA11
  float	d = dot(E0, D);                   // fB0
  float	e = dot(E1, D);                   // fB1
  float	f = dot(D, D);                    // fC
  float det = fabsf(a*c - b*b);
  s = b*e-c*d;
  t = b*d-a*e;
	
  OpenMesh::Vec3f     residual;

  float distance2;

  if ( s + t <= det )
  {
    if ( s < 0.0f )
    {
      if ( t < 0.0f )  // region 4
      {
        if ( d < 0.0f )
        {
          t = 0.0f;
          if ( -d >= a )
          {
            s = 1.0f;
            distance2 = a+2.0f*d+f;
          }
          else
          {
            s = -d/a;
            distance2 = d*s+f;
          }
        }
        else
        {
          s = 0.0f;
          if ( e >= 0.0f )
          {
            t = 0.0f;
            distance2 = f;
          }
          else if ( -e >= c )
          {
            t = 1.0f;
            distance2 = c+2.0f*e+f;
          }
          else
          {
            t = -e/c;
            distance2 = e*t+f;
          }
        }
      }
      else  // region 3
      {
        s = 0.0f;
        if ( e >= 0.0f )
        {
          t = 0.0f;
          distance2 = f;
        }
        else if ( -e >= c )
        {
          t = 1.0f;
          distance2 = c+2.0f*e+f;
        }
        else
        {
          t = -e/c;
          distance2 = e*t+f;
        }
      }
    }
    else if ( t < 0.0f )  // region 5
    {
      t = 0.0f;
      if ( d >= 0.0f )
      {
        s = 0.0f;
        distance2 = f;
      }
      else if ( -d >= a )
      {
        s = 1.0f;
        distance2 = a+2.0f*d+f;
      }
      else
      {
        s = -d/a;
        distance2 = d*s+f;
      }
    }
    else  // region 0
    {
      // minimum at interior point
      float inv_det = 1.0f/det;
      s *= inv_det;
      t *= inv_det;
      distance2 = s*(a*s+b*t+2.0f*d) +
        t*(b*s+c*t+2.0f*e)+f;
    }
  }
  else
  {
    float tmp0, tmp1, numer, denom;

    if ( s < 0.0f )  // region 2
    {
      tmp0 = b + d;
      tmp1 = c + e;
      if ( tmp1 > tmp0 )
      {
        numer = tmp1 - tmp0;
        denom = a-2.0f*b+c;
        if ( numer >= denom )
        {
          s = 1.0f;
          t = 0.0f;
          distance2 = a+2.0f*d+f;
        }
        else
        {
          s = numer/denom;
          t = 1.0f - s;
          distance2 = s*(a*s+b*t+2.0f*d) +
            t*(b*s+c*t+2.0f*e)+f;
        }
      }
      else
      {
        s = 0.0f;
        if ( tmp1 <= 0.0f )
        {
          t = 1.0f;
          distance2 = c+2.0f*e+f;
        }
        else if ( e >= 0.0f )
        {
          t = 0.0f;
          distance2 = f;
        }
        else
        {
          t = -e/c;
          distance2 = e*t+f;
        }
      }
    }
    else if ( t < 0.0f )  // region 6
    {
      tmp0 = b + e;
      tmp1 = a + d;
      if ( tmp1 > tmp0 )
      {
        numer = tmp1 - tmp0;
        denom = a-2.0f*b+c;
        if ( numer >= denom )
        {
          t = 1.0f;
          s = 0.0f;
          distance2 = c+2.0f*e+f;
        }
        else
        {
          t = numer/denom;
          s = 1.0f - t;
          distance2 = s*(a*s+b*t+2.0f*d)+ t*(b*s+c*t+2.0f*e)+f;
        }
      }
      else
      {
        t = 0.0f;
        if ( tmp1 <= 0.0f )
        {
          s = 1.0f;
          distance2 = a+2.0f*d+f;
        }
        else if ( d >= 0.0f )
        {
          s = 0.0f;
          distance2 = f;
        }
        else
        {
          s = -d/a;
          distance2 = d*s+f;
        }
      }
    }
    else  // region 1
    {
      numer = c + e - b - d;
      if ( numer <= 0.0f )
      {
        s = 0.0f;
        t = 1.0f;
        distance2 = c+2.0f*e+f;
      }
      else
      {
        denom = a-2.0f*b+c;
        if ( numer >= denom )
        {
          s = 1.0f;
          t = 0.0f;
          distance2 = a+2.0f*d+f;
        }
        else
        {
          s = numer/denom;
          t = 1.0f - s;
          distance2 = s*(a*s+b*t+2.0f*d) + t*(b*s+c*t+2.0f*e)+f;
        }
      }
    }
  }

  residual = p - (B + s*E0 + t*E1);

  return	residual;
}

// ============================================================================