size_t closestPointToRay(const V3f* points, size_t nPoints,
                         const V3f& rayOrigin, const V3f& rayDirection,
                         double longitudinalScale, double* distance)
{
    const V3f T = rayDirection.normalized();
    const double f = longitudinalScale*longitudinalScale - 1;
    size_t nearestIdx = -1;
    double nearestDist2 = DBL_MAX;
    for(size_t i = 0; i < nPoints; ++i)
    {
        const V3f v = points[i] - rayOrigin; // vector from ray origin to point
        const double a = v.dot(T); // distance along ray to point of closest approach to test point
        const double r2 = v.length2() + f*a*a;

        if(r2 < nearestDist2)
        {
            // new closest angle to axis
            nearestDist2 = r2;
            nearestIdx = i;
        }
    }
    if(distance)
    {
        if(nPoints == 0)
            *distance = DBL_MAX;
        else
            *distance = sqrt(nearestDist2);
    }
    return nearestIdx;
}

	bool intersectP(CRay& _ray) {
		if (_ray.m_id == m_id)
			return false; 
		V3f oc = _ray.m_pos - m_c; 
		float a = _ray.m_dir.dot(_ray.m_dir); // dir should be unit
		float b = _ray.m_dir.dot(oc);
		float c = oc.dot(oc) - m_r2; 
		float delta = b * b  - a * c; 
		if (delta < 0)   // no solution 
			return false; 
		else if (delta > -EPS && delta < EPS) {  // one solution
			float t = - b / a;
			if (t > _ray.m_t_max || t < _ray.m_t_min)  // out of range
				return false; 
		} else {   // two solutions 
			float deltasqrt = sqrt(delta);
			float t1 = (- b - deltasqrt) / a;
			float t2 = (- b + deltasqrt) / a;
			if (t1 >= _ray.m_t_max || t1 <= _ray.m_t_min)
				return false; 

			if (t2 >= _ray.m_t_max || t2 <= _ray.m_t_min)
				return false; 
		}

		return true; 
	}

void PhongBrdf::randVonMisesFisher3(V3f mu, float kappa, int n, V3f* directions) {


	V3f normal(0,0,1);
	V3f u = mu.cross(normal);
	float cost = dot(mu,normal);
	float sint = u.length();
	u = u.normalize();

	M33f rot(cost + u.x * u.x * (1 - cost),
			u.x * u.y * (1 - cost) - u.z * sint,
			u.x * u.z * (1 - cost) + u.y * sint,
			u.y * u.x * (1 - cost) + u.z * sint,
			cost + u.y * u.y * (1 - cost),
			u.y * u.z * (1 - cost) - u.x * sint,
			u.z * u.x * (1 - cost) - u.y * sint,
			u.z	* u.y * (1 - cost) + u.x * sint,
			cost + u.z * u.z * (1 - cost));

	float c = 2/kappa*(sinh(kappa)); // normalizing constant

	float y, w, v;
	for (int i=0; i < n; i++) {
		y = randomGenerator.RandomFloat();
		w = 1/kappa * log( exp(-kappa) + kappa * c * y );
		v = 2*M_PI*randomGenerator.RandomFloat();

		directions[i].x = sqrt(1-w*w)*cos(v);
		directions[i].y = sqrt(1-w*w)*sin(v);
		directions[i].z = w;

		directions[i] = directions[i]*rot;
	}

}

void PhongBrdf::getRandomDirections(const V3f incoming,
			const V3f normal, int n, V3f* directions) {

	V3f R = -incoming - 2 * (dot(-incoming, normal)) * normal;
	R = R.normalize();

	randVonMisesFisher3(R, phongExponent, n, directions);
}

//-*****************************************************************************
void MeshDrwHelper::updateNormals( V3fArraySamplePtr iN )
{
    if ( !m_valid || !m_meshP )
    {
        makeInvalid();
        return;
    }

//std::cout << "normals - " << m_name << std::endl;

    // Now see if we need to calculate normals.
    if ( ( m_meshN && iN == m_meshN ) )//||
//         ( !iN && m_customN.size() > 0 ) )
    {
        return;
    }

    size_t numPoints = m_meshP->size();
    m_meshN = iN;
    m_customN.clear();

    // Right now we only handle "vertex varying" normals,
    // which have the same cardinality as the points
    if ( !m_meshN || m_meshN->size() != numPoints )
    {
        // Make some custom normals.
        m_meshN.reset();
        m_customN.resize( numPoints );
        std::fill( m_customN.begin(), m_customN.end(), V3f( 0.0f ) );

        //std::cout << "Recalcing normals for object: "
        //          << m_host.name() << std::endl;

        for ( size_t tidx = 0; tidx < m_triangles.size(); ++tidx )
        {
            const Tri &tri = m_triangles[tidx];

            const V3f &A = (*m_meshP)[tri[0]];
            const V3f &B = (*m_meshP)[tri[1]];
            const V3f &C = (*m_meshP)[tri[2]];

            V3f AB = B - A;
            V3f AC = C - A;

            V3f wN = AB.cross( AC );
            m_customN[tri[0]] += wN;
            m_customN[tri[1]] += wN;
            m_customN[tri[2]] += wN;
        }

        // Normalize normals.
        for ( size_t nidx = 0; nidx < numPoints; ++nidx )
        {
            m_customN[nidx].normalize();
        }
    }
}

//Splice edge with centered sphere.
bool SpliceEdgeWithSphere(const V3f & a, const V3f & b, float radius, V3f * out){

  float al = a.length(); 
  float bl = b.length();
  
  if( ( (al >= radius) && (bl >= radius) ) ||
      ( (al <= radius) && (bl <= radius) ) ) return false; 

  *out =  a * (bl-radius)/(bl-al) + b * (radius-al)/(bl-al);

  return true;
};

void
Camera::rotateVectorQuat(float angle, float x, float y, float z, V3f& vector) const
{
	bm::quaternion<float> temp(x * sin(angle / 2.0f), y * sin(angle / 2.0f), z
			* sin(angle / 2.0f), cos(angle / 2.0f));
	bm::quaternion<float> quat_vec(*vector.x, *vector.y, *vector.z, 0.0);
	bm::quaternion<float> result = (temp * quat_vec) * bm::conj(temp);

	result.real();
	vector.setX(result.R_component_1());
	vector.setY(result.R_component_2());
	vector.setZ(result.R_component_3());
}

void
Camera::mouseRotate(float angleY, float angleZ)
{
	V3f vAxis = V3f::cross(mView - mEye, mUp);
	vAxis.normalize();

	// Rotate around our perpendicular axis and along the y-axis
//	rotateView(angleZ, vAxis);
//	rotateView(angleY, 0.0f, 1.0f, 0.0f);
	rotateView(angleZ, vAxis);
	rotateView(angleY, 0.0f, 1.0f,0.0f);
//	applyToGL();
}

V2f	
latLong (const V3f &dir)
{
    float r = sqrt (dir.z * dir.z + dir.x * dir.x);

    float latitude = (r < abs (dir.y))?
			 acos (r / dir.length()) * sign (dir.y):
			 asin (dir.y / dir.length());

    float longitude = (dir.z == 0 && dir.x == 0)? 0: atan2 (dir.x, dir.z);

    return V2f (latitude, longitude);
}

int main(){
  V3f x(0,0,1); V3f xr(rot_x(x, 0.87)); same("x rotation", x.dot(xr), cos(0.87));
  V3f y(0,0,1); V3f yr(rot_y(y, 0.23)); same("y rotation", y.dot(yr), cos(0.23));
  V3f z(1,0,0); V3f zr(rot_z(z, 0.19)); same("z rotation", z.dot(zr), cos(0.19));

  V3f nx(3,2,5);
  V3f ny(-2,3,4);
  V3f nz(-4,4,3.8);

  V3f nnx(3,2,5);
  V3f nny(-2,3,4);
  V3f nnz(-4,4,3.8);

  ortoNormalize(nnx, nny, nnz);
  
  same("x unit", nnx.length(), 1.0);
  same("y unit", nny.length(), 1.0);
  same("z unit", nnz.length(), 1.0);

  V3f tmp; tmp.cross(nnx, nx);

  same("x colinear", tmp.length(), 0.0);
  
  tmp.cross(nnx, nny); tmp-=nnz; same("x orto", tmp.length(), 0);
  tmp.cross(nny, nnz); tmp-=nnx; same("y orto", tmp.length(), 0);
  tmp.cross(nnz, nnx); tmp-=nny; same("z orto", tmp.length(), 0);


};

bool SpherePrimitiveEvaluator::closestPoint( const V3f &p, PrimitiveEvaluator::Result *result ) const
{
	assert( dynamic_cast<Result *>( result ) );

	Result *sr = static_cast<Result *>( result );

	sr->m_p = p.normalized() * m_sphere->radius();

	return true;
}

//make basis ortonormal again.
void ortoNormalize(V3f & nnx, V3f & nny, V3f & nnz){
  V3f newy; newy.cross(nnz, nnx);
  V3f newz; newz.cross(nnx, newy);
  newy /= newy.length();
  newz /= newz.length();
  nnx /= nnx.length();
  nny = newy;
  nnz = newz;
};

	bool intersect(CRay& _ray, float* _thit, CLocalGeo* _local, int& _id)  {
		if (_ray.m_id == m_id)
			return false; 
		float t; 
		V3f oc = _ray.m_pos - m_c; 
		float a = _ray.m_dir.dot(_ray.m_dir); // dir should be unit
		float b = _ray.m_dir.dot(oc);
		float c = oc.dot(oc) - m_r2; 
		float delta = b * b  - a * c; 
		if (delta < 0)   // no solution 
			return false; 
		else if (delta > -EPS && delta < EPS) {  // one solution
			t = - b / a;
			if (t > _ray.m_t_max || t < _ray.m_t_min)  // out of range
				return false; 
		} else {   // two solutions 
			float deltasqrt = sqrt(delta);
			float t1 = (- b - deltasqrt) / a;
			float t2 = (- b + deltasqrt) / a;
			bool flag = false; 
			t = _ray.m_t_max; 
			if (t1 <= _ray.m_t_max && t1 >= _ray.m_t_min) {
				flag = true; 
				t = min(t, t1);
			}

			if (t2 <= _ray.m_t_max && t2 >= _ray.m_t_min) {
				flag = true; 
				t = min(t, t2);
			}
			if (!flag)   // both out of range
				return false; 
		}

		// pass t, compute CLocalGeo
		*_thit = t; 
		_id = m_id; 
		_local->m_pos = _ray.Ray_t(t);
		_local->m_n = _local->m_pos - m_c; 
		_local->m_n = _local->m_n / _local->m_n.norm(); 
		return true; 
	} 

void PhongModelApprox::approximate(const Hemisphere& hemi) {

	N = hemi.getNormal();
	phong = hemi.getPhong();
	V3f* directions = hemi.getLobeDirections();
	C3f* radiosities = hemi.getLobeRadiosities();


	for (int i=0; i <lobeDirs.size(); i++) {
		if (i >= hemi.getNLobes()) {
			float nan = std::numeric_limits<float>::quiet_NaN();
			lobeDirs[i] = V3f(nan,nan,nan);
		} else {
			V3f L = directions[i];
			L = L/L.length();
			V3f R = -L - 2 * (dot(-L, N)) * N;

			lobeDirs[i] = R;
			lobeCols[i] = radiosities[i]*dot(N,L);
		}
	}
}

void CameraController::tumble( const Imath::V2f &p )
{
	V2f d = p - m_data->motionStart;

	V3f centreOfInterestInWorld = V3f( 0, 0, -m_data->centreOfInterest ) * m_data->motionMatrix;
	V3f xAxisInWorld = V3f( 1, 0, 0 );
	m_data->motionMatrix.multDirMatrix( xAxisInWorld, xAxisInWorld );
	xAxisInWorld.normalize();

	M44f t;
	t.translate( centreOfInterestInWorld );

		t.rotate( V3f( 0, -d.x / 100.0f, 0 ) );

		M44f xRotate;
		xRotate.setAxisAngle( xAxisInWorld, -d.y / 100.0f );

		t = xRotate * t;

	t.translate( -centreOfInterestInWorld );

	m_data->transform->matrix = m_data->motionMatrix * t;
}

int ImagePrimitiveEvaluator::intersectionPoints( const V3f &origin, const V3f &direction,
                std::vector<PrimitiveEvaluator::ResultPtr> &results, float maxDistance ) const
{
	results.clear();

	V3f hitPoint;
	Box3f bound = m_image->bound();
	bool hit = boxIntersects( bound , origin, direction.normalized(), hitPoint );

	if ( hit )
	{
		if ( ( origin - hitPoint ).length2() < maxDistance * maxDistance )
		{
			ResultPtr result = staticPointerCast< Result >( createResult() );
			result->m_p = hitPoint;

			results.push_back( result );
		}
	}

	return results.size();
}

static void renderDiskExact(IntegratorT& integrator, V3f p, V3f n, float r) {

    int faceRes = integrator.res();
    float plen2 = p.length2();
    if(plen2 == 0) // Sanity check
        return;
    // Angle from face normal to edge is acos(1/sqrt(3)).
    static float cosFaceAngle = 1.0f/sqrtf(3);
    static float sinFaceAngle = sqrtf(2.0f/3.0f);

    // iterate over all the faces.
    for(int iface = MicroBuf::Face_begin; iface < MicroBuf::Face_begin; ++iface) {

        // Cast this back to a Face enum.
        MicroBuf::Face face = static_cast<MicroBuf::Face>(iface);

        // Avoid rendering to the current face if the disk definitely doesn't
        // touch it.  First check the cone angle
        if(sphereOutsideCone(p, plen2, r, MicroBuf::faceNormal(face),
                             cosFaceAngle, sinFaceAngle))
            continue;
        float dot_pFaceN = MicroBuf::dotFaceNormal(face, p);
        float dot_nFaceN = MicroBuf::dotFaceNormal(face, n);
        // If the disk is behind the camera and the disk normal is relatively
        // aligned with the face normal (to within the face cone angle), the
        // disk can't contribute to the face and may be culled.
        if(dot_pFaceN < 0 && fabs(dot_nFaceN) > cosFaceAngle)
            continue;
        // Check whether disk spans the perspective divide for the current
        // face.  Note: sin^2(angle(n, faceN)) = (1 - dot_nFaceN*dot_nFaceN)
        if((1 - dot_nFaceN*dot_nFaceN)*r*r >= dot_pFaceN*dot_pFaceN)
        {
            // When the disk spans the perspective divide, the shape of the
            // disk projected onto the face is a hyperbola.  Bounding a
            // hyperbola is a pain, so the easiest thing to do is trace a ray
            // for every pixel on the face and check whether it hits the disk.
            //
            // Note that all of the tricky rasterization rubbish further down
            // could probably be replaced by the following ray tracing code if
            // I knew a way to compute the tight raster bound.
            integrator.setFace(face);
            for(int iv = 0; iv < faceRes; ++iv)
                for(int iu = 0; iu < faceRes; ++iu)
                {
                    // V = ray through the pixel
                    V3f V = integrator.rayDirection(face, iu, iv);
                    // Signed distance to plane containing disk
                    float t = dot(p, n)/dot(V, n);
                    if(t > 0 && (t*V - p).length2() < r*r)
                    {
                        // The ray hit the disk, record the hit
                        integrator.addSample(iu, iv, t, 1.0f);
                    }
                }
            continue;
        }
        // If the disk didn't span the perspective divide and is behind the
        // camera, it may be culled.
        if(dot_pFaceN < 0)
            continue;
        // Having gone through all the checks above, we know that the disk
        // doesn't span the perspective divide, and that it is in front of the
        // camera.  Therefore, the disk projected onto the current face is an
        // ellipse, and we may compute a quadratic function
        //
        //   q(u,v) = a0*u*u + b0*u*v + c0*v*v + d0*u + e0*v + f0
        //
        // such that the disk lies in the region satisfying q(u,v) < 0.  To do
        // this, start with the implicit definition of the disk on the plane,
        //
        //   norm(dot(p,n)/dot(V,n) * V - p)^2 - r^2 < 0
        //
        // and compute coefficients A,B,C such that
        //
        //   A*dot(V,V) + B*dot(V,n)*dot(p,V) + C < 0
        float dot_pn = dot(p,n);
        float A = dot_pn*dot_pn;
        float B = -2*dot_pn;
        float C = plen2 - r*r;
        // Project onto the current face to compute the coefficients a0 through
        // to f0 for q(u,v)
        V3f pp = MicroBuf::canonicalFaceCoords(face, p);
        V3f nn = MicroBuf::canonicalFaceCoords(face, n);
        float a0 = A + B*nn.x*pp.x + C*nn.x*nn.x;
        float b0 = B*(nn.x*pp.y + nn.y*pp.x) + 2*C*nn.x*nn.y;
        float c0 = A + B*nn.y*pp.y + C*nn.y*nn.y;
        float d0 = (B*(nn.x*pp.z + nn.z*pp.x) + 2*C*nn.x*nn.z);
        float e0 = (B*(nn.y*pp.z + nn.z*pp.y) + 2*C*nn.y*nn.z);
        float f0 = (A + B*nn.z*pp.z + C*nn.z*nn.z);
        // Finally, transform the coefficients so that they define the
        // quadratic function in *raster* face coordinates, (iu, iv)
        float scale = 2.0f/faceRes;
        float scale2 = scale*scale;
        float off = 0.5f*scale - 1.0f;
        float a = scale2*a0;
        float b = scale2*b0;
        float c = scale2*c0;
        float d = ((2*a0 + b0)*off + d0)*scale;
        float e = ((2*c0 + b0)*off + e0)*scale;
        float f = (a0 + b0 + c0)*off*off + (d0 + e0)*off + f0;
//.........这里部分代码省略.........

void
ObjModelLoader::secondPass()
{
	string lastType = "";
	int grpId = -1;
	mPriModelPtr->setCurrentGrp(mPriModelPtr->getGrpStart()->first);

	fs::path my_path(mPriFName);
	fs::ifstream objFile;
	objFile.open(my_path, ios::out);
	char line[200];
	vector<string> tokens;
	while (objFile.getline(line, 200)) {
		tokens = splitSpace(string(line));
		if (tokens[0] == ("g")) {
//			removeSpecialCharsFromName(tokens[1]);
			++grpId;
			lastType = "g";
			string longname("");
			for (unsigned int i = 1; i < tokens.size(); ++i) {
				longname.append(tokens[i]);
				longname.append("_");
			}
			longname.erase(longname.end() - 1);

			mPriModelPtr->setCurrentGrp(longname);
			//cout << _model.getCurrentGrpPtr()->name << endl;
		}
		else if (tokens[0] == ("v")) {
			lastType = "v";
		}
		// moved material assignment to 2nd pass because it references a group.
		// But only at end of 1st pass we know our groups
		else if (tokens[0] == ("usemtl")) {
			//			cout << "found material reference " << tokens[1] << " in obj-file" << endl;
			string longName("");
			for (unsigned int i = 1; i < tokens.size(); ++i) {
				longName.append(tokens[i]);
				longName.append(" ");
			}
			longName.erase(longName.end() - 1);
			if (mPriMatMap.find(longName) != mPriMatMap.end())
				mPriModelPtr->getCurrentGrpPtr()->setMat(*mPriMatMap[longName]);
			lastType = "usemtl";
		}
		else if (tokens[0] == ("f")) {
			//cout << "Number of Components per face: " << tokens.size()-1 << endl;
			Face* f = new Face();
			f->norm = false;
			f->vert = false;
			f->matIdx = 0;
			f->tex = 0;
			f->fNormal = 0;

			string::size_type loc = tokens[1].find("/", 0);
			if (loc != string::npos) {
				for (string::size_type i = 1; i < tokens.size(); ++i) {
					vector<int> comp = extractFaceComponents(tokens[i]);
					// vertices
					if ((comp[0] & 4)) {
						V3f* vtx = mPriModelPtr->getVPtr(comp[1] - 1);
						f->vertexPtrList.push_back(vtx);
						vtx->addFaceRef(f);
						f->vert = true;
						mPriModelPtr->getCurrentGrpPtr()->nVertices++;
						mPriModelPtr->incVCount();
						mPriModelPtr->getCurrentGrpPtr()->bb->expand(*vtx);
					}
					// textures
					if ((comp[0] & 2)) {
						f->texturePtrList.push_back(mPriModelPtr->getTPtr(
								comp[2] - 1));
						++f->tex;
						mPriModelPtr->getCurrentGrpPtr()->nTextureCoords++;
						mPriModelPtr->incTCount();
					}
					// normals
					if ((comp[0] & 1)) {
						f->normalPtrList.push_back(mPriModelPtr->getNPtr(
								comp[3] - 1));
						f->norm = true;
						mPriModelPtr->getCurrentGrpPtr()->nNormals++;
						mPriModelPtr->incNCount();
					}
					comp.clear();
				}
			}
			else {
				V3f* vtx = mPriModelPtr->getVPtr(atoi(tokens[1].c_str()) - 1);
				f->vertexPtrList.push_back(vtx);
				mPriModelPtr->getCurrentGrpPtr()->bb->expand(*vtx);
				vtx = mPriModelPtr->getVPtr(atoi(tokens[2].c_str()) - 1);
				f->vertexPtrList.push_back(vtx);
				mPriModelPtr->getCurrentGrpPtr()->bb->expand(*vtx);
				vtx = mPriModelPtr->getVPtr(atoi(tokens[3].c_str()) - 1);
				f->vertexPtrList.push_back(vtx);
				mPriModelPtr->getCurrentGrpPtr()->bb->expand(*vtx);
				f->vert = true;
				mPriModelPtr->getCurrentGrpPtr()->nVertices += 3;
				mPriModelPtr->incVCount(3);
//.........这里部分代码省略.........

//-*****************************************************************************
void MeshDrwHelper::draw( const DrawContext & iCtx ) const
{
    // Bail if invalid.
    if ( !m_valid || m_triangles.size() < 1 || !m_meshP )
    {
        return;
    }

    const V3f *points = m_meshP->get();
    const V3f *normals = NULL;
    if ( m_meshN  && ( m_meshN->size() == m_meshP->size() ) )
    {
        normals = m_meshN->get();
    }
    else if ( m_customN.size() == m_meshP->size() )
    {
        normals = &(m_customN.front());
    }

#ifndef SIMPLE_ABC_VIEWER_NO_GL_CLIENT_STATE
//#if 0
    {
        GL_NOISY( glEnableClientState( GL_VERTEX_ARRAY ) );
        if ( normals )
        {
            GL_NOISY( glEnableClientState( GL_NORMAL_ARRAY ) );
            GL_NOISY( glNormalPointer( GL_FLOAT, 0,
                                       ( const GLvoid * )normals ) );
        }

        GL_NOISY( glVertexPointer( 3, GL_FLOAT, 0,
                                   ( const GLvoid * )points ) );

        GL_NOISY( glDrawElements( GL_TRIANGLES,
                                  ( GLsizei )m_triangles.size() * 3,
                                  GL_UNSIGNED_INT,
                                  ( const GLvoid * )&(m_triangles[0]) ) );

        if ( normals )
        {
            GL_NOISY( glDisableClientState( GL_NORMAL_ARRAY ) );
        }
        GL_NOISY( glDisableClientState( GL_VERTEX_ARRAY ) );
    }
#else
    glBegin( GL_TRIANGLES );

    for ( size_t i = 0; i < m_triangles.size(); ++i )
    {
        const Tri &tri = m_triangles[i];
        const V3f &vertA = points[tri[0]];
        const V3f &vertB = points[tri[1]];
        const V3f &vertC = points[tri[2]];

        if ( normals )
        {
            const V3f &normA = normals[tri[0]];
            glNormal3fv( ( const GLfloat * )&normA );
            glVertex3fv( ( const GLfloat * )&vertA );

            const V3f &normB = normals[tri[1]];
            glNormal3fv( ( const GLfloat * )&normB );
            glVertex3fv( ( const GLfloat * )&vertB );

            const V3f &normC = normals[tri[2]];
            glNormal3fv( ( const GLfloat * )&normC );
            glVertex3fv( ( const GLfloat * )&vertC );
        }
        else
        {
            V3f AB = vertB - vertA;
            V3f AC = vertC - vertA;
            V3f N = AB.cross( AC );
            if ( N.length() > 1.0e-4f )
            {
                N.normalize();
                glNormal3fv( ( const GLfloat * )&N );
            }

            glVertex3fv( ( const GLfloat * )&vertA );

            glVertex3fv( ( const GLfloat * )&vertB );

            glVertex3fv( ( const GLfloat * )&vertC );
        }

    }

    glEnd();

#endif
}

std::pair<PrimitiveVariable, PrimitiveVariable> IECoreScene::MeshAlgo::calculateTangents(
	const MeshPrimitive *mesh,
	const std::string &uvSet, /* = "uv" */
	bool orthoTangents, /* = true */
	const std::string &position /* = "P" */
)
{
	if( mesh->minVerticesPerFace() != 3 || mesh->maxVerticesPerFace() != 3 )
	{
		throw InvalidArgumentException( "MeshAlgo::calculateTangents : MeshPrimitive must only contain triangles" );
	}

	const V3fVectorData *positionData = mesh->variableData<V3fVectorData>( position );
	if( !positionData )
	{
		std::string e = boost::str( boost::format( "MeshAlgo::calculateTangents : MeshPrimitive has no Vertex \"%s\" primitive variable." ) % position );
		throw InvalidArgumentException( e );
	}

	const V3fVectorData::ValueType &points = positionData->readable();

	const IntVectorData *vertsPerFaceData = mesh->verticesPerFace();
	const IntVectorData::ValueType &vertsPerFace = vertsPerFaceData->readable();

	const IntVectorData *vertIdsData = mesh->vertexIds();
	const IntVectorData::ValueType &vertIds = vertIdsData->readable();

	const auto uvIt = mesh->variables.find( uvSet );
	if( uvIt == mesh->variables.end() || uvIt->second.interpolation != PrimitiveVariable::FaceVarying || uvIt->second.data->typeId() != V2fVectorDataTypeId )
	{
		throw InvalidArgumentException( ( boost::format( "MeshAlgo::calculateTangents : MeshPrimitive has no FaceVarying V2fVectorData primitive variable named \"%s\"."  ) % ( uvSet ) ).str() );
	}

	const V2fVectorData *uvData = runTimeCast<V2fVectorData>( uvIt->second.data.get() );
	const V2fVectorData::ValueType &uvs = uvData->readable();

	// I'm a little unsure about using the vertIds as a fallback for the stIndices.
	const IntVectorData::ValueType &uvIndices = uvIt->second.indices ? uvIt->second.indices->readable() : vertIds;

	size_t numUVs = uvs.size();

	std::vector<V3f> uTangents( numUVs, V3f( 0 ) );
	std::vector<V3f> vTangents( numUVs, V3f( 0 ) );
	std::vector<V3f> normals( numUVs, V3f( 0 ) );

	for( size_t faceIndex = 0; faceIndex < vertsPerFace.size(); faceIndex++ )
	{
		assert( vertsPerFace[faceIndex] == 3 );

		// indices into the facevarying data for this face
		size_t fvi0 = faceIndex * 3;
		size_t fvi1 = fvi0 + 1;
		size_t fvi2 = fvi1 + 1;
		assert( fvi2 < vertIds.size() );
		assert( fvi2 < uvIndices.size() );

		// positions for each vertex of this face
		const V3f &p0 = points[vertIds[fvi0]];
		const V3f &p1 = points[vertIds[fvi1]];
		const V3f &p2 = points[vertIds[fvi2]];

		// uv coordinates for each vertex of this face
		const V2f &uv0 = uvs[uvIndices[fvi0]];
		const V2f &uv1 = uvs[uvIndices[fvi1]];
		const V2f &uv2 = uvs[uvIndices[fvi2]];

		// compute tangents and normal for this face
		const V3f e0 = p1 - p0;
		const V3f e1 = p2 - p0;

		const V2f e0uv = uv1 - uv0;
		const V2f e1uv = uv2 - uv0;

		V3f tangent = ( e0 * -e1uv.y + e1 * e0uv.y ).normalized();
		V3f bitangent = ( e0 * -e1uv.x + e1 * e0uv.x ).normalized();

		V3f normal = ( p2 - p1 ).cross( p0 - p1 );
		normal.normalize();

		// and accumlate them into the computation so far
		uTangents[uvIndices[fvi0]] += tangent;
		uTangents[uvIndices[fvi1]] += tangent;
		uTangents[uvIndices[fvi2]] += tangent;

		vTangents[uvIndices[fvi0]] += bitangent;
		vTangents[uvIndices[fvi1]] += bitangent;
		vTangents[uvIndices[fvi2]] += bitangent;

		normals[uvIndices[fvi0]] += normal;
		normals[uvIndices[fvi1]] += normal;
		normals[uvIndices[fvi2]] += normal;

	}

	// normalize and orthogonalize everything
	for( size_t i = 0; i < uTangents.size(); i++ )
	{
		normals[i].normalize();

		uTangents[i].normalize();
//.........这里部分代码省略.........