C++ Vec3r::norm方法代码示例

void Grid::castInfiniteRay(const Vec3r& orig,
			   const Vec3r& dir,
			   OccludersSet& occluders,
			   unsigned timestamp) {
  Vec3r end = Vec3r(orig + FLT_MAX * dir / dir.norm());
  bool inter = initInfiniteRay(orig, dir, timestamp);
  if(!inter)
      return;
  allOccludersGridVisitor visitor(occluders);
  castRayInternal(visitor);
}

inline static real angle(WOEdge *h)
{
	const Vec3r& n1 = h->GetbFace()->GetNormal();
	const Vec3r& n2 = h->GetaFace()->GetNormal();
	const Vec3r v = h->GetVec();
	real sine = (n1 ^ n2) * v / v.norm();
	if (sine >= 1.0) {
		return M_PI / 2.0;
	}
	if (sine <= -1.0) {
		return -M_PI / 2.0;
	}
	return ::asin(sine);
}

/*! gts_vertex_principal_directions:
 *  @v: a #WVertex.
 *  @s: a #GtsSurface.
 *  @Kh: mean curvature normal (a #Vec3r).
 *  @Kg: Gaussian curvature (a real).
 *  @e1: first principal curvature direction (direction of largest curvature).
 *  @e2: second principal curvature direction.
 *
 *  Computes the principal curvature directions at a point given @Kh and @Kg, the mean curvature normal and
 *  Gaussian curvatures at that point, computed with gts_vertex_mean_curvature_normal() and
 *  gts_vertex_gaussian_curvature(), respectively.
 *
 *  Note that this computation is very approximate and tends to be unstable. Smoothing of the surface or the principal
 *  directions may be necessary to achieve reasonable results.
 */
void gts_vertex_principal_directions(WVertex *v, Vec3r Kh, real Kg, Vec3r &e1, Vec3r &e2)
{
	Vec3r N;
	real normKh;

	Vec3r basis1, basis2, d, eig;
	real ve2, vdotN;
	real aterm_da, bterm_da, cterm_da, const_da;
	real aterm_db, bterm_db, cterm_db, const_db;
	real a, b, c;
	real K1, K2;
	real *weights, *kappas, *d1s, *d2s;
	int edge_count;
	real err_e1, err_e2;
	int e;
	WVertex::incoming_edge_iterator itE;

	/* compute unit normal */
	normKh = Kh.norm();

	if (normKh > 0.0) {
		Kh.normalize();
	}
	else {
		/* This vertex is a point of zero mean curvature (flat or saddle point). Compute a normal by averaging
		 * the adjacent triangles
		 */
		N[0] = N[1] = N[2] = 0.0;

		for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++)
			N = Vec3r(N + (*itE)->GetaFace()->GetNormal());
		real normN = N.norm();
		if (normN <= 0.0)
			return;
		N.normalize();
	}

	/* construct a basis from N: */
	/* set basis1 to any component not the largest of N */
	basis1[0] =  basis1[1] =  basis1[2] = 0.0;
	if (fabs (N[0]) > fabs (N[1]))
		basis1[1] = 1.0;
	else
		basis1[0] = 1.0;

	/* make basis2 orthogonal to N */
	basis2 = (N ^ basis1);
	basis2.normalize();

	/* make basis1 orthogonal to N and basis2 */
	basis1 = (N ^ basis2);
	basis1.normalize();

	aterm_da = bterm_da = cterm_da = const_da = 0.0;
	aterm_db = bterm_db = cterm_db = const_db = 0.0;
	int nb_edges = v->GetEdges().size();

	weights = (real *)malloc(sizeof(real) * nb_edges);
	kappas = (real *)malloc(sizeof(real) * nb_edges);
	d1s = (real *)malloc(sizeof(real) * nb_edges);
	d2s = (real *)malloc(sizeof(real) * nb_edges);
	edge_count = 0;

	for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) {
		WOEdge *e;
		WFace *f1, *f2;
		real weight, kappa, d1, d2;
		Vec3r vec_edge;
		if (!*itE)
			continue;
		e = *itE;

		/* since this vertex passed the tests in gts_vertex_mean_curvature_normal(), this should be true. */
		//g_assert(gts_edge_face_number (e, s) == 2);

		/* identify the two triangles bordering e in s */
		f1 = e->GetaFace();
		f2 = e->GetbFace();

		/* We are solving for the values of the curvature tensor
		 *     B = [ a b ; b c ].
		 *  The computations here are from section 5 of [Meyer et al 2002].
		 *
		 *  The first step is to calculate the linear equations governing the values of (a,b,c). These can be computed
		 *  by setting the derivatives of the error E to zero (section 5.3).
//.........这里部分代码省略.........

bool FEdge::intersectParametric(FEdge & fe2, Vec3r viewpoint, real t3D, real u3D)
{
  Vec3r A1 = vertexA()->getPoint3D();
  Vec3r B1 = vertexB()->getPoint3D();
  Vec3r A2 = fe2.vertexA()->getPoint3D();
  Vec3r B2 = fe2.vertexB()->getPoint3D();

  if (sameSide(A1,B1,viewpoint, A2, B2) || sameSide(A2, B2, viewpoint, A1, B1))
    return false;

  // now, there *must* be an intersection.

  // for each edge, the normal of the plane containing the edge and the viewpoint
  Vec3r N1 = (A1-viewpoint) ^ (B1-viewpoint);
  Vec3r N2 = (A2-viewpoint) ^ (B2-viewpoint);

  // direction vector of the intersection of the two planes.
  Vec3r V = N1 ^ N2;

  // check if the planes coincide (i.e., source edges are colinear)
  assert(V.norm() > 0);
  
  // ----- compute t parameter ------

  // form a plane for line 1, normal to the plane containing the viewpoint

  Vec3r BA1 = B1 - A1;
  Vec3r hsNormal1 = N1 ^ BA1;
  
  // intersect ray in direction of V through the plane
  real w1;
  GeomUtils::intersection_test res1 = GeomUtils::intersectLinePlanePN(viewpoint, V, hsNormal1, A1, w1);

  if (res1 != GeomUtils::DO_INTERSECT)
    {
      printf("res1 = %d\n", res1);
      printf("viewpoint = [%f %f %f]\n", viewpoint[0], viewpoint[1], viewpoint[2]);
      printf("A1 = [%f %f %f]\n", A1[0], A1[1], A1[2]);
      printf("B1 = [%f %f %f]\n", B1[0], B1[1], B1[2]);
      printf("A2 = [%f %f %f]\n", A2[0], A2[1], A2[2]);
      printf("B2 = [%f %f %f]\n", B2[0], B2[1], B2[2]);
      printf("N1 = [%f %f %f]\n", N1[0], N1[1], N1[2]);
      printf("N2 = [%f %f %f]\n", N2[0], N2[1], N2[2]);
      printf("V = [%f %f %f]\n", V[0], V[1], V[2]);
      printf("hsNormal1 = [%f %f %f]\n", hsNormal1[0], hsNormal1[1], hsNormal1[2]);
    }

  assert(res1 == GeomUtils::DO_INTERSECT);

  Vec3r pt1 = viewpoint + w1 * V;

  t3D = ((pt1 - A1) * BA1) / (BA1*BA1);

  assert(t3D >=0 && t3D <= 1);

  // if (t3D < 0 || t3D > 1)
    //    return false;


  // ----- compute u parameter ------

  // form a half-space plane for line 2

  Vec3r BA2 = B2 - A2;
  Vec3r hsNormal2 = N2 ^ BA2;
  
  real w2;
  GeomUtils::intersection_test res2 = GeomUtils::intersectLinePlanePN(viewpoint, V, hsNormal2, A2, w2);

  if (res2 != GeomUtils::DO_INTERSECT)
    {
      printf("res1 = %d\n", res1);
      printf("viewpoint = [%f %f %f]\n", viewpoint[0], viewpoint[1], viewpoint[2]);
      printf("A1 = [%f %f %f]\n", A1[0], A1[1], A1[2]);
      printf("B1 = [%f %f %f]\n", B1[0], B1[1], B1[2]);
      printf("A2 = [%f %f %f]\n", A2[0], A2[1], A2[2]);
      printf("B2 = [%f %f %f]\n", B2[0], B2[1], B2[2]);
      printf("N1 = [%f %f %f]\n", N1[0], N1[1], N1[2]);
      printf("N2 = [%f %f %f]\n", N2[0], N2[1], N2[2]);
      printf("V = [%f %f %f]\n", V[0], V[1], V[2]);
      printf("hsNormal2 = [%f %f %f]\n", hsNormal2[0], hsNormal2[1], hsNormal2[2]);
    }

  assert(res2 == GeomUtils::DO_INTERSECT);

  Vec3r pt2 = viewpoint + w2 * V;

  u3D = ((pt2 - A2) * BA2) / (BA2*BA2);

  assert( u3D >=0 && u3D <=1);

  //  if (u3D < 0 || u3D > 1)
  //    return false;


  
  return true;
}