C++ vec3::cross方法代码示例

 VRC(const vec3& VRP, const vec3& VPN, const vec3& VUP)
 {
   O = VRP;
   n = VPN;
   u = VUP.cross(n).versor();
   v = VPN.cross(u);
 }

    static void lookAt(const vec3<Type>& eye, const vec3<Type>& center, const vec3<Type>& up)
    {
        const vec3<Type> forward = (center - eye).normalized();
        const vec3<Type> side = forward.cross(up).normalized();
        const vec3<Type> upVector = side.cross(forward);

        matrix4<Type> m;
        m.m_data[0][0] = side.x();
        m.m_data[1][0] = side.y();
        m.m_data[2][0] = side.z();
        m.m_data[3][0] = -side.dot(eye);
        m.m_data[0][1] = upVector.x();
        m.m_data[1][1] = upVector.y();
        m.m_data[2][1] = upVector.z();
        m.m_data[3][1] = -upVector.dot(eye);
        m.m_data[0][2] = -forward.x();
        m.m_data[1][2] = -forward.y();
        m.m_data[2][2] = -forward.z();
        m.m_data[3][2] = -forward.dot(eye);
        m.m_data[0][3] = 0;
        m.m_data[1][3] = 0;
        m.m_data[2][3] = 0;
        m.m_data[3][3] = 1;

        return m;
    }

void	OGLViewMatrix(const vec3& dir, const vec3& up, const vec3& loc)
// Multiplies the current OpenGL matrix by the view matrix corresponding to the given parameters.
{
	matrix m;
	vec3 Z = dir.cross(up);
	m.SetColumn(0, -Z);
	m.SetColumn(1, up);
	m.SetColumn(2, dir);
	m.SetColumn(3, loc);
	m.Invert();
	float	mat[16];

	// Copy to the 4x4 layout.
	for (int col = 0; col < 4; col++) {
		for (int row = 0; row < 3; row++) {
			mat[col * 4 + row] = m.GetColumn(col).Get(row);
		}
		if (col < 3) {
			mat[col * 4 + 3] = 0;
		} else {
			mat[col * 4 + 3] = 1;
		}
	}

	// Apply to the current OpenGL matrix.
	glMultMatrixf(mat);
}

void
RayTracer::renderImage(Image& image, bool isAdaptative)
//[]---------------------------------------------------[]
//|  Run the ray tracer                                 |
//[]---------------------------------------------------[]
{
	clock_t t = clock();

	image.getSize(W, H);
	// init auxiliary VRC
	VRC_n = camera->getViewPlaneNormal();
	VRC_v = camera->getViewUp();
	VRC_u = VRC_v.cross(VRC_n);
	// init auxiliary mapping variables
	I_w = Math::inverse<REAL>(REAL(W));
	I_h = Math::inverse<REAL>(REAL(H));

	REAL height = camera->windowHeight();

	W >= H ? V_w = (V_h = height) * W * I_h : V_h = (V_w = height) * H * I_w;
	if (isAdaptative)
		adaptativeScan(image);
	else
		scan(image);
	printf("\nNumber of rays: %lu", numberOfRays);
	printf("\nNumber of hits: %lu", numberOfHits);
	printElapsedTime("\nDONE! ", clock() - t);
}

  const std::pair<vec3, vec3> rotate(const vec3 &r, const vec3 &vt, const real dt)
  {
    const real th = vt.abs()*dt/r.abs();

    const real sinth = sin(th);
    const real costh = cos(th);

    const vec3 rv = r.cross(vt);
    const vec3 n = rv/rv.abs();
    const vec3 n2(sqr(n.x), sqr(n.y), sqr(n.z));

    const real Axx = n2.x + (n2.y + n2.z)*costh;
    const real Ayy = n2.y + (n2.x + n2.z)*costh;
    const real Azz = n2.z + (n2.x + n2.y)*costh;
    const real Axy = n.x*n.y*(1-costh) - n.z*sinth;
    const real Axz = n.x*n.z*(1-costh) + n.y*sinth; 
    const real Ayx = n.x*n.y*(1-costh) + n.z*sinth;
    const real Ayz = n.y*n.z*(1-costh) - n.x*sinth;
    const real Azx = n.x*n.z*(1-costh) - n.y*sinth;
    const real Azy = n.y*n.z*(1-costh) + n.x*sinth;

    return std::make_pair(
        vec3(
          vec3(Axx, Axy, Axz)*r,
          vec3(Ayx, Ayy, Ayz)*r,
          vec3(Azx, Azy, Azz)*r),
        vec3(
          vec3(Axx, Axy, Axz)*vt,
          vec3(Ayx, Ayy, Ayz)*vt,
          vec3(Azx, Azy, Azz)*vt)
        );
  }

void camera::MoveTo(vec3 p0, vec3 pref, vec3 up) {
  vec3 n = (p0 - pref).normalize();
  vec3 u = up.cross(n).normalize();
  vec3 v = n.cross(u);

  mat4 M;
  M[0]  = u.x;
  M[1]  = u.y;
  M[2]  = u.z;
  M[4]  = v.x;
  M[5]  = v.y;
  M[6]  = v.z;
  M[8]  = n.x;
  M[9]  = n.y;
  M[10] = n.z;
  M[15] = 1;
  mat4 M_inv = M.Transpose();

  mat4 T      = mat4::Identity().Translate(vec3(0, 0, 0) - p0);
  mat4 T_inv  = mat4::Identity().Translate(p0);

  // Create new view and create inverse view for 
  //  later delta calculations
  view          = M * T;
  view_inverse  = T_inv * M_inv;

  // Cache camera for parametric curving
  cached_p0 = p0;
  cached_pref = pref;
  cached_up = up;
}

void Camera::setSide(const vec3& s)
{
    vec3 u = up();
    vec3 d = normalize(s.cross(u));
    vec3 p = position();
    vec3 e(-dot(s, p), -dot(u, p), dot(d, p));
    _modelViewMatrix[0] = vec4(s.x, u.x, -d.x, 0.0f);
    _modelViewMatrix[1] = vec4(s.y, u.y, -d.y, 0.0f);
    _modelViewMatrix[2] = vec4(s.z, u.z, -d.z, 0.0f);
    _modelViewMatrix[3] = vec4(e.x, e.y,  e.z, 1.0f);
    modelViewUpdated();
}

void mat4::lookAt(const vec3 &eye, const vec3 &center, const vec3 &up)
{
	vec3 z = vec3(eye-center).getNormal();
	vec3 y = up;
	vec3 x = up.cross(z);

	identity();
	setAxisX(x);
	setAxisY(y);
	setAxisZ(z);
	setPos(eye);
}

void	mouse_motion_handler(float dx, float dy, int state)
// Handles mouse motion.
{
	bool	left_button = (state & 1) ? true : false;
	bool	right_button = (state & 4) ? true : false;

	float	f_speed = (1 << speed) / 200.f;

	if (left_button && right_button) {
		// Translate in the view plane.
		vec3	right = viewer_dir.cross(viewer_up);
		viewer_pos += right * dx * f_speed + viewer_up * -dy * f_speed;
		
	} else if (left_button) {
		// Rotate the viewer.
		viewer_theta += -dx / 100;
		while (viewer_theta < 0) viewer_theta += (float) (2 * M_PI);
		while (viewer_theta >= 2 * M_PI) viewer_theta -= (float) (2 * M_PI);

		viewer_phi += -dy / 100;
		const float	plimit = (float) (M_PI / 2);
		if (viewer_phi > plimit) viewer_phi = plimit;
		if (viewer_phi < -plimit) viewer_phi = -plimit;

		viewer_dir = vec3(1, 0, 0);
		viewer_up = vec3(0, 1, 0);
		
		viewer_dir = Geometry::Rotate(viewer_theta, viewer_up, viewer_dir);
		vec3	right = viewer_dir.cross(viewer_up);
		viewer_dir = Geometry::Rotate(viewer_phi, right, viewer_dir);
		viewer_up = Geometry::Rotate(viewer_phi, right, viewer_up);

	} else if (right_button) {
		// Translate the viewer.
		vec3	right = viewer_dir.cross(viewer_up);
		viewer_pos += right * dx * f_speed + viewer_dir * -dy * f_speed;
	}
}

/**
 * @brief LookAt Setup view matrix
 * @param eye Camera position
 * @param center Looking center
 * @param up Up vector
 * @return
 */
mat4 LookAt(vec3 eye, vec3 view, vec3 up)
{
    vec3 z_cam = view.normalized();
    vec3 x_cam = up.cross(z_cam).normalized();
    vec3 y_cam = z_cam.cross(x_cam);
    mat3 R = mat3::Zero();
    R.col(0) = x_cam;
    R.col(1) = y_cam;
    R.col(2) = z_cam;
    mat4 L = mat4::Zero();
    L.block(0, 0, 3, 3) =  R.transpose();
    L(3, 3) = 1.0f;
    L.block(0, 3, 3, 1) = -R.transpose() * (eye);
    return L;
}

        /*! Intersect the ray with the given triangle defined by vert0,vert1,vert2.
        Return true if there is an intersection.
        If there is an intersection, a vector a_tuv is returned, where t is the
        distance to the plane in which the triangle lies and (u,v) represents the
        coordinates inside the triangle.
        */
        bool intersect(const vec3<Type>& vert0, const vec3<Type>& vert1, const vec3<Type>& vert2, vec3<Type>& a_tuv) const
        {
            // Tomas Moller and Ben Trumbore.
            // Fast, minimum storage ray-triangle intersection.
            // Journal of graphics tools, 2(1):21-28, 1997

            // find vectors for two edges sharing vert0
            const vec3<Type> edge1 = vert1 - vert0;
            const vec3<Type> edge2 = vert2 - vert0;

            // begin calculating determinant - also used to calculate U parameter
            const vec3<Type> pvec = m_direction.cross(edge2);

            // if determinant is near zero, ray lies in plane of triangle
            const Type det = edge1.dot(pvec);

            if (det < EPS)
                return false;

            // calculate distance from vert0 to ray origin
            const vec3<Type> tvec = m_origin - vert0;

            // calculate U parameter and test bounds
            a_tuv[1] = tvec.dot(pvec);

            if (a_tuv[1] < 0.0 || a_tuv[1] > det)
                return false;

            // prepare to test V parameter
            const vec3<Type> qvec = tvec.cross(edge1);

            // calculate V parameter and test bounds
            a_tuv[2] = m_origin.dot(qvec);

            if (a_tuv[2] < 0.0 || a_tuv[1] + a_tuv[2] > det)
                return false;

            // calculate t, scale parameters, ray intersects triangle
            a_tuv[0] = edge2.dot(qvec);

            const Type inv_det = (Type)1.0 / det;

            a_tuv[0] *= inv_det;
            a_tuv[1] *= inv_det;
            a_tuv[2] *= inv_det;

            return true;
        }

void Shadow::calculateFrustumVertices(const Light &light, const Actor &actor, float lx, float ly, vec3 &ntl, vec3 &ntr, vec3 &nbl, vec3 &nbr, vec3 &ftl, vec3 &ftr, vec3 &fbl, vec3 &fbr)
{
	const vec3 &lightPosition = light.getPosition();
	const vec3 &lightCenter = actor.getPos();
	const vec3 lightUp(0.0f, 1.0f, 0.0f);

	float tanX = tanf(lx);
	float tanY = tanf(ly);
	float nh = nearClip * tanY;
	float nw = nearClip * tanX;
	float fh = farClip  * tanY;
	float fw = farClip  * tanX;

	vec3 dir,nc,fc,X,Y,Z;

	Z = lightPosition - lightCenter;
	Z.normalize();

	X = lightUp.cross(Z);
	X.normalize();

	Y = Z.cross(X);

	nc = lightPosition - Z * nearClip;
	fc = lightPosition - Z * farClip;

	ntl = nc + Y * nh - X * nw;
	ntr = nc + Y * nh + X * nw;
	nbl = nc - Y * nh - X * nw;
	nbr = nc - Y * nh + X * nw;

	ftl = fc + Y * fh - X * fw;
	ftr = fc + Y * fh + X * fw;
	fbl = fc - Y * fh - X * fw;
	fbr = fc - Y * fh + X * fw;
}

	void system::compute_fluid_grad(const bool do_ngb)
	{
#if 0
		const int nimport = Uimport.size(); 
		assert(nimport == (int)site_import.size());
		Wextra_import.resize(nimport);

    const int nactive_site = site_active_list.size();
		const int nactive_site_with_ngb = nactive_site + (do_ngb ? site_with_ngb_active_list.size() : 0);
		
    for (int isite = 0; isite < nactive_site_with_ngb; isite++)
		{
			const int i = isite < nactive_site ? 
				site_active_list[isite] : site_with_ngb_active_list[isite - nactive_site];

			const Cell     &ci     = cell_list[i];
			const Fluid_st &Wst_i  = W_st  [i];
			const Fluid    &Wi     = Wst_i.w;
			const vec3     &ipos   = Wst_i.pos;

			vec3 Ji(0.0);     // current,  J = Del x B

			const vec3 Bi(Wi[Fluid::BX], Wi[Fluid::BY], Wi[Fluid::BZ]);
			
			const int nface = ci.faces().size();	
			for (int iface = 0; iface < nface; iface++)
			{
				const Face &face = face_list[ci.faces()[iface]];

				vec3 dri = face.centroid - ipos;
        if      (dri.x >  0.5*global_domain_size.x) dri.x -= global_domain_size.x;
        else if (dri.x < -0.5*global_domain_size.x) dri.x += global_domain_size.x;
        if      (dri.y >  0.5*global_domain_size.y) dri.y -= global_domain_size.y;
        else if (dri.y < -0.5*global_domain_size.y) dri.y += global_domain_size.y;
        if      (dri.z >  0.5*global_domain_size.z) dri.z -= global_domain_size.z;
        else if (dri.z < -0.5*global_domain_size.z) dri.z += global_domain_size.z;
        assert(std::abs(dri.x) < 0.5*global_domain_size.x);
        assert(std::abs(dri.y) < 0.5*global_domain_size.y);
        assert(std::abs(dri.z) < 0.5*global_domain_size.z);
        
        const vec3 centroid = dri;
				const real area     = face.area();
				const vec3 normal   = face.n * ((centroid * face.n < 0.0) ? (-1.0/area) : (1.0/area));

				const real dsh  = centroid * normal;
				assert(dsh > 0.0);

				const real ids  = (dsh != 0.0) ? 0.5/dsh : 0.0;
				const real idsA = area * ids;

				const vec3 drh = normal * (dsh * idsA);
				const vec3 fij = (centroid - drh) * idsA;
				const int j = site_map(face.ngb<false>(i));
				assert(j >= 0);

				const Fluid &Wj = W_st[j].w;
        const vec3 Bj(Wj[Fluid::BX], Wj[Fluid::BY], Wj[Fluid::BZ]);
        Ji += drh.cross(Bj+Bi) + fij.cross(Bj-Bi);
			}

			const real invV = 1.0/cell_list[i].Volume;
			Ji *= ptcl_import[i].etaOhm * invV;
      
      if (do_ngb) W_st[i].J = Ji;
      else     
      {
        const real dt = t_global - std::abs(W_st[i].tlast);
//        assert(dt > 0.0);
//        W_rec[i].J = (Ji - W_st[i].J) * (1.0/dt);
      }

		}
#endif	
	}

inline vec3
vAxis(const vec3& DOP, const vec3& VUP)
{
  vec3 u = DOP.cross(VUP).versor();
  return u.cross(DOP);
}

void EditorToolBar::onLeftMouseDown()
{
	ASSERT(world!=0, "world was null!  Call setWorld first!");

	if(g_GUI.mouseOverSomeWidget) return;

	// Get a position on the ground beneath the cursor
	float e = selected ? selected->getPos().y : 0.0f;
	const vec3 groundPos = getGroundPickPos(e);

	// Grab the pool of objects we are working from
	ActorSet &objects = world->getObjects();

	// Get the id of the object under the mouse cursor
	OBJECT_ID id = objects.getClosest<Actor>(groundPos, 2.0f);

	// Process the action depending on the current tool
	switch(toolBarTools->getTool())
	{
	case ToolBarForEditorTools::EDITOR_SELECT_TOOL:
		if(objects.isMember(id))
		{
			Actor * p = &objects.get(id);
			showActorPane(p);
		}
		else
		{
			hideActorPane();
		}
		break;

	case ToolBarForEditorTools::EDITOR_MOVE_TOOL:
		{
			if(selected)
			{
				// Hold the right mouse button to slide objects along the y-axis
				if(g_Input.MouseRight)
				{
					vec3 delta = groundPos - selected->getPos();
					delta.y=0;
					float dist = delta.getMagnitude()*0.3f;

					// Place the object on this spot
					selected->Place(vec3(selected->getPos().x, dist, selected->getPos().z));
				}
				else
				{
					selected->Place(groundPos);
				}
			}
		}
		break;

	case ToolBarForEditorTools::EDITOR_ROTATE_TOOL:
		{
			if(selected)
			{
				const vec3 delta = vec3(selected->getPos().x-groundPos.x, 0, selected->getPos().z-groundPos.z);
				const vec3 zAxis = delta.getNormal();
				const vec3 yAxis = vec3(0,1,0);
				const vec3 xAxis = yAxis.cross(zAxis).getNormal();

				mat4 orientation = selected->getOrientation();
				orientation.setAxisZ(zAxis);
				orientation.setAxisY(yAxis);
				orientation.setAxisX(xAxis);
				orientation.setPos(vec3(0,0,0));
				selected->setOrientation(orientation);
			}
		}
		break;

	case ToolBarForEditorTools::EDITOR_ROTATE_X_TOOL:
		{
			if(selected)
			{
				vec3 delta = groundPos - selected->getPos();
				float angle = atan2f(delta.x, delta.y) * 0.1f;

				mat4 rot;
				rot.rotateX(angle);

		                mat4 orientation = selected->getOrientation();

				orientation *= rot;

				selected->setOrientation(orientation);
			}
		}
		break;

	case ToolBarForEditorTools::EDITOR_ROTATE_Z_TOOL:
		{
			if(selected)
			{
				vec3 delta = groundPos - selected->getPos();
				float angle = atan2f(delta.x, delta.y) * 0.1f;

				mat4 rot;
				rot.rotateZ(angle);
//.........这里部分代码省略.........