C++ Vector_3::squared_length方法代码示例

static std::vector<Point_3> mapPointsFromOXYplane(std::vector<Point_2> points,
		Vector_3 nu)
{
	DEBUG_START;

	ASSERT(!!Vector3d(nu.x(), nu.y(), nu.z()) && "nu is null vector");

	Vector_3 ez(0., 0, 1.);
	double length = sqrt(nu.squared_length());
	ASSERT(std::fpclassify(length) != FP_ZERO);
	nu = nu * 1. / length; /* Normalize std::vector \nu. */
	ASSERT(std::isfinite(nu.x()));
	ASSERT(std::isfinite(nu.y()));
	ASSERT(std::isfinite(nu.z()));
	Vector_3 tau = cross_product(nu, ez);

	std::vector<Point_3> pointsMapped;
	CGAL::Origin o;
	for (auto &point : points)
	{
		pointsMapped.push_back(o + tau * point.x() + ez * point.y());
	}
	DEBUG_END;
	return pointsMapped;
}

double geometryUtils::computeVoronoiArea(Vertex_handle vertex) {
    double voronoiArea = 0.0;
    Vertex_circulator j;

    j = vertex->vertex_begin();

    do {
        Point_3 p11 = j->vertex()->point();
        Point_3 p12 = j->next()->vertex()->point();
        Point_3 p13 = j->next()->next()->vertex()->point();
        Vector_3 v11 = p13 - p12;
        Vector_3 v12 = p11 - p12;
        v11 = v11 / sqrt(CGAL::to_double(v11.squared_length()));
        v12 = v12 / sqrt(CGAL::to_double(v12.squared_length()));
        double alpha = acos(CGAL::to_double(v11 * v12));
        Point_3 p22 = j->opposite()->vertex()->point();
        Point_3 p23 = j->opposite()->next()->vertex()->point();
        Vector_3 v21 = p11 - p23;
        Vector_3 v22 = p22 - p23;
        v21 = v21 / sqrt(CGAL::to_double(v21.squared_length()));
        v22 = v22 / sqrt(CGAL::to_double(v22.squared_length()));

        double beta = acos(CGAL::to_double(v21 * v22));
        Vector_3 x = p13 - p11;
        double length = CGAL::to_double(x.squared_length());

        voronoiArea += (1.0 / 8.0) * (1.0 / tan(alpha) + 1.0 / tan(beta)) * length;

    } while (++j != vertex->vertex_begin());

    return voronoiArea;

};

void add_neighbourhood_to_hullPoints(pointVector& hullPoints, const Vertex_const_handle& vert, const double& tapeSize) {
	Point_3 pnt = vert->point();
	hullPoints.push_back(pnt);								// Add vertex point
	Nef_polyhedron::SVertex_const_iterator svcIt = vert->svertices_begin(), svcItEND = vert->svertices_end();
	CGAL_For_all(svcIt,svcItEND) {
		Vector_3 vecR(pnt,svcIt->target()->point());
		Vector_3 vecRnew = vecR * tapeSize / std::sqrt(CGAL::to_double(vecR.squared_length()));
		if ((vecR.squared_length()-OVERLAP_DIST_THRESHOLD) > vecRnew.squared_length())
			hullPoints.push_back(pnt+vecRnew);						// Add svertex neighbourhood point (tapesize away from vertex)
		else
			hullPoints.push_back(svcIt->target()->point());
	}

void buildFacets(Polyhedron_3 polyhedron,
		std::vector<SimpleEdge_3> &edges,
		std::vector<Vector_3> &U, std::vector<double> &H,
		std::map<int, int> &indices)
{
	DEBUG_START;
	std::vector<Plane_3> planes = renumerateFacets(polyhedron, edges,
			indices);
	for (const Plane_3 &plane : planes)
	{
		Vector_3 norm = plane.orthogonal_vector();
		double length = sqrt(norm.squared_length());
		norm = norm * (1. / length);
		double value = -plane.d() / length;
		ASSERT(value > 0.);
		U.push_back(norm);
		H.push_back(value);
	}
	DEBUG_END;
}

void buildMainTopology(std::vector<SimpleEdge_3> &edges,
		std::vector<Vector_3> &u,
		std::vector<double> &h, std::vector<Vector_3> &points,
		FixedTopology *FT)
{
	DEBUG_START;
	unsigned iTangient = 0;
	for (unsigned i = 0; i < edges.size(); ++i)
	{
		Vector_3 A = edges[i].A;
		Vector_3 B = edges[i].B;
		for (const Plane_3 plane : edges[i].tangients)
		{
			Vector_3 norm = plane.orthogonal_vector();
			double length = sqrt(norm.squared_length());
			norm = norm * (1. / length);
			double value = -plane.d() / length;
			ASSERT(value > 0.);
			u.push_back(norm);
			h.push_back(value);

			if (norm * A > norm * B)
				FT->tangient[2 * i].insert(iTangient);
			else
				FT->tangient[2 * i + 1].insert(iTangient);
			++iTangient;
		}

		FT->incident[edges[i].iForward].insert(2 * i);
		FT->incident[edges[i].iForward].insert(2 * i + 1);
		FT->incident[edges[i].iBackward].insert(2 * i);
		FT->incident[edges[i].iBackward].insert(2 * i + 1);
		points.push_back(A);
		points.push_back(B);
	}
	ASSERT(iTangient > 0);
	DEBUG_END;
}

void Viewer::drawEdge(const Point_3& from, const Point_3& to, const QColor& clr, float r)
{
  /* Draw regular lines */
  if( m_isFlat ) {
    // disable lighting
    ::glDisable( GL_LIGHTING );

    ::glLineWidth(1.0);
  	qglColor( clr );

    ::glBegin(GL_LINES);
    ::glVertex3f( from.x(), from.y(), from.z() );
    ::glVertex3f( to.x(), to.y(), to.z() );
    ::glEnd();

    // resume lighting
    ::glEnable( GL_LIGHTING );

    return;
  }

  /* Draw edges as 3D cylinders */
  GLboolean lighting, colorMaterial;
  ::glGetBooleanv( GL_LIGHTING, &lighting );
  ::glGetBooleanv( GL_COLOR_MATERIAL, &colorMaterial );
  ::glEnable( GL_LIGHTING );
  ::glDisable(GL_COLOR_MATERIAL);

  float color[4];
  color[0] = clr.redF();
  color[1] = clr.greenF();
  color[2] = clr.blueF();
  color[3] = clr.alphaF();

  Vector_3 v = to - from;

  // compute the length of the edge
  // method 1:
//  float length = sqrt( CGAL::squared_distance( from, to ) );
  // method 2:
  float length = sqrt( v.squared_length() );

  // normalize
  v = v / length;
  // compute the angle: cos theta = v.z/1.0
  GLfloat angle = acos( v.z() ) / 3.1415927 * 180;

  ::glPushMatrix();

  // move to "from" point
  ::glTranslatef( from.x(), from.y(), from.z() );
  // rotate from z-axis to from-->to
  //  axis: cross product of z-axis and from-->to
  ::glRotatef( angle, -v.y(), v.x(), 0.0f );
  // draw
  GLUquadricObj* quadratic = ::gluNewQuadric();	// Create A Pointer To The Quadric Object
  ::gluQuadricNormals( quadratic, GLU_SMOOTH );	// Create Smooth Normals
  ::glMaterialfv( GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, color );
  // gluCylinder draws a cylinder oriented along the z-axis
  ::gluCylinder( quadratic, r, r, length, 16, 4 );

  // move back to origin
  ::glPopMatrix();

  if ( colorMaterial )
    ::glEnable( GL_COLOR_MATERIAL );
  if ( !lighting )
    ::glDisable( GL_LIGHTING );
}

//.........这里部分代码省略.........
  {
    assert(!S1.is_degenerate());
    // All points on S2 are the same
    P2 = S2.source();
    point_segment_squared_distance(P2,S1,P1,dst);
    return true;
  }

  assert(!S1.is_degenerate());
  assert(!S2.is_degenerate());

  Vector_3 u = S1.target() - S1.source();
  Vector_3 v = S2.target() - S2.source();
  Vector_3 w = S1.source() - S2.source();

  const auto a = u.dot(u);         // always >= 0
  const auto b = u.dot(v);
  const auto c = v.dot(v);         // always >= 0
  const auto d = u.dot(w);
  const auto e = v.dot(w);
  const auto D = a*c - b*b;        // always >= 0
  assert(D>=0);
  const auto sc=D, sN=D, sD = D;       // sc = sN / sD, default sD = D >= 0
  const auto tc=D, tN=D, tD = D;       // tc = tN / tD, default tD = D >= 0

  bool parallel = false;
  // compute the line parameters of the two closest points
  if (D==0) 
  { 
    // the lines are almost parallel
    parallel = true;
    sN = 0.0;         // force using source point on segment S1
    sD = 1.0;         // to prevent possible division by 0.0 later
    tN = e;
    tD = c;
  } else
  {
    // get the closest points on the infinite lines
    sN = (b*e - c*d);
    tN = (a*e - b*d);
    if (sN < 0.0) 
    { 
      // sc < 0 => the s=0 edge is visible
      sN = 0.0;
      tN = e;
      tD = c;
    } else if (sN > sD) 
    {  // sc > 1  => the s=1 edge is visible
      sN = sD;
      tN = e + b;
      tD = c;
    }
  }

  if (tN < 0.0) 
  {
    // tc < 0 => the t=0 edge is visible
    tN = 0.0;
    // recompute sc for this edge
    if (-d < 0.0)
    {
      sN = 0.0;
    }else if (-d > a)
    {
      sN = sD;
    }else 
    {
      sN = -d;
      sD = a;
    }
  }else if (tN > tD) 
  {
    // tc > 1  => the t=1 edge is visible
    tN = tD;
    // recompute sc for this edge
    if ((-d + b) < 0.0)
    {
      sN = 0;
    }else if ((-d + b) > a)
    {
      sN = sD;
    }else
    {
      sN = (-d +  b);
      sD = a;
    }
  }
  // finally do the division to get sc and tc
  sc = (abs(sN) == 0 ? 0.0 : sN / sD);
  tc = (abs(tN) == 0 ? 0.0 : tN / tD);

  // get the difference of the two closest points
  P1 = S1.source() + sc*(S1.target()-S1.source());
  P2 = S2.source() + sc*(S2.target()-S2.source());
  Vector_3   dP = w + (sc * u) - (tc * v);  // =  S1(sc) - S2(tc)
  assert(dP == (P1-P2));

  dst = dP.squared_length();   // return the closest distance
  return parallel;
}