Initial commit:

PhysX 3.4.0 Update @ 21294896
APEX 1.4.0 Update @ 21275617

[CL 21300167]

1 files changed, 1006 insertions, 0 deletions

diff --git a/APEX_1.4/module/clothing/embedded/ExtClothGeodesicTetherCooker.cpp b/APEX_1.4/module/clothing/embedded/ExtClothGeodesicTetherCooker.cpp
new file mode 100644
index 00000000..e0fea857
--- /dev/null
+++ b/APEX_1.4/module/clothing/embedded/ExtClothGeodesicTetherCooker.cpp
@@ -0,0 +1,1006 @@
+/*
+ * Copyright (c) 2008-2015, NVIDIA CORPORATION.  All rights reserved.
+ *
+ * NVIDIA CORPORATION and its licensors retain all intellectual property
+ * and proprietary rights in and to this software, related documentation
+ * and any modifications thereto.  Any use, reproduction, disclosure or
+ * distribution of this software and related documentation without an express
+ * license agreement from NVIDIA CORPORATION is strictly prohibited.
+ */
+
+// Copyright (c) 2004-2008 AGEIA Technologies, Inc. All rights reserved.
+// Copyright (c) 2001-2004 NovodeX AG. All rights reserved.  
+
+#include "ExtClothConfig.h"
+#if APEX_USE_CLOTH_API
+
+#include "ExtClothTetherCooker.h"
+#include "PxStrideIterator.h"
+
+// from shared foundation
+
+#include <PsArray.h>
+#include <PsSort.h>
+#include <Ps.h>
+#include <PsMathUtils.h>
+#include "PxVec4.h"
+#include "PsIntrinsics.h"
+
+using namespace nvidia;
+using namespace physx;
+
+namespace
+{
+	// calculate the inclusive prefix sum, equivalent of std::partial_sum
+	template <typename T>
+	void prefixSum(const T* first, const T* last, T* dest)
+	{
+		if (first != last)
+		{	
+			*(dest++) = *(first++);
+			for (; first != last; ++first, ++dest)
+				*dest = *(dest-1) + *first;
+		}
+	}
+
+	template <typename T>
+	void gatherAdjacencies(shdfnd::Array<uint32_t>& valency, shdfnd::Array<uint32_t>& adjacencies, 
+		const PxBoundedData& triangles, const PxBoundedData& quads)
+	{
+		// count number of edges per vertex
+		PxStrideIterator<const T> tIt, qIt;
+		tIt = physx::PxMakeIterator((const T*)triangles.data, triangles.stride);
+		for(uint32_t i=0; i<triangles.count; ++i, ++tIt, ++qIt)
+		{
+			for(uint32_t j=0; j<3; ++j)
+				valency[tIt.ptr()[j]] += 2;
+		}
+		qIt = physx::PxMakeIterator((const T*)quads.data, quads.stride);
+		for(uint32_t i=0; i<quads.count; ++i, ++tIt, ++qIt)
+		{
+			for(uint32_t j=0; j<4; ++j)
+				valency[qIt.ptr()[j]] += 2;
+		}
+
+		prefixSum(valency.begin(), valency.end(), valency.begin());
+		adjacencies.resize(valency.back());
+
+		// gather adjacent vertices
+		tIt = physx::PxMakeIterator((const T*)triangles.data, triangles.stride);
+		for(uint32_t i=0; i<triangles.count; ++i, ++tIt)
+		{
+			for(uint32_t j=0; j<3; ++j)
+			{
+				adjacencies[--valency[tIt.ptr()[j]]] = tIt.ptr()[(j+1)%3];
+				adjacencies[--valency[tIt.ptr()[j]]] = tIt.ptr()[(j+2)%3];
+			}
+		}
+		qIt = physx::PxMakeIterator((const T*)quads.data, quads.stride);
+		for(uint32_t i=0; i<quads.count; ++i, ++qIt)
+		{
+			for(uint32_t j=0; j<4; ++j)
+			{
+				adjacencies[--valency[qIt.ptr()[j]]] = qIt.ptr()[(j+1)%4];
+				adjacencies[--valency[qIt.ptr()[j]]] = qIt.ptr()[(j+3)%4];
+			}
+		}
+	}
+
+	template <typename T>
+	void gatherIndices(shdfnd::Array<uint32_t>& indices, 
+		const PxBoundedData& triangles, const PxBoundedData& quads)
+	{
+		PxStrideIterator<const T> tIt, qIt;
+
+		indices.reserve(triangles.count * 3 + quads.count * 6);
+	
+		tIt = physx::PxMakeIterator((const T*)triangles.data, triangles.stride);
+		for(uint32_t i=0; i<triangles.count; ++i, ++tIt)
+		{
+			indices.pushBack(tIt.ptr()[0]);
+			indices.pushBack(tIt.ptr()[1]);
+			indices.pushBack(tIt.ptr()[2]);
+		}
+		qIt = physx::PxMakeIterator((const T*)quads.data, quads.stride);
+		for(uint32_t i=0; i<quads.count; ++i, ++qIt)
+		{
+			indices.pushBack(qIt.ptr()[0]);
+			indices.pushBack(qIt.ptr()[1]);
+			indices.pushBack(qIt.ptr()[2]);
+			indices.pushBack(qIt.ptr()[0]);
+			indices.pushBack(qIt.ptr()[2]);
+			indices.pushBack(qIt.ptr()[3]);
+		}
+	}
+
+	// maintain heap status after elements have been pushed (heapify)
+	template<typename T>
+	void pushHeap(shdfnd::Array<T> &heap, const T &value)
+	{
+		heap.pushBack(value);
+		T* begin = heap.begin();
+		T* end = heap.end();
+
+		if (end <= begin)
+			return;
+	
+		uint32_t current = uint32_t(end - begin) - 1;
+		while (current > 0)
+		{
+			const uint32_t parent = (current - 1) / 2;
+			if (!(begin[parent] < begin[current]))
+				break;
+
+			shdfnd::swap(begin[parent], begin[current]);
+			current = parent;
+		}
+	}
+
+	// pop one element from the heap
+	template<typename T>
+	T popHeap(shdfnd::Array<T> &heap)
+	{
+		T* begin = heap.begin();
+		T* end = heap.end();
+
+		shdfnd::swap(begin[0], end[-1]); // exchange elements
+
+		// shift down
+		end--;
+
+		uint32_t current = 0;
+		while (begin + (current * 2 + 1) < end)
+		{
+			uint32_t child = current * 2 + 1;
+			if (begin + child + 1 < end && begin[child] < begin[child + 1])
+				++child;
+
+			if (!(begin[current] < begin[child]))
+				break;
+
+			shdfnd::swap(begin[current], begin[child]);
+			current = child;
+		}
+
+		return heap.popBack();
+	}
+
+	// ---------------------------------------------------------------------------------------
+	struct VertexDistanceCount
+	{
+		VertexDistanceCount(int vert, float dist, int count) 
+			: vertNr(vert), distance(dist), edgeCount(count) {}
+
+		int vertNr;
+		float distance;
+		int edgeCount;
+		bool operator < (const VertexDistanceCount& v) const
+		{
+			return v.distance < distance;
+		}
+	};
+
+	// ---------------------------------------------------------------------------------------
+	struct PathIntersection
+	{
+		uint32_t vertOrTriangle;
+		uint32_t index; // vertex id or triangle edge id
+		float s; // only used for edge intersection
+		float distance; // computed distance
+
+	public:
+		PathIntersection() {}
+
+		PathIntersection(uint32_t vort, uint32_t in_index, float in_distance, float in_s = 0.0f) 
+			: vertOrTriangle(vort), index(in_index), s(in_s), distance(in_distance)
+		{
+		}
+	};
+
+	//---------------------------------------------------------------------------------------
+	struct VertTriangle
+	{
+		VertTriangle(int vert, int triangle)
+			: mVertIndex(vert), mTriangleIndex(triangle)
+		{
+		}
+
+		bool operator<(const VertTriangle &vt) const
+		{
+			return mVertIndex == vt.mVertIndex ?
+				mTriangleIndex < vt.mTriangleIndex : mVertIndex < vt.mVertIndex;
+		}
+
+		int mVertIndex;
+		int mTriangleIndex;
+	};
+
+	// ---------------------------------------------------------------------------------------
+	struct MeshEdge
+	{
+		MeshEdge(int v0, int v1, int halfEdgeIndex)
+			: mFromVertIndex(v0), mToVertIndex(v1), mHalfEdgeIndex(halfEdgeIndex)
+		{
+			if(mFromVertIndex > mToVertIndex)
+				shdfnd::swap(mFromVertIndex, mToVertIndex);
+		}
+
+		bool operator<(const MeshEdge& e) const
+		{
+			return mFromVertIndex == e.mFromVertIndex ?
+				mToVertIndex < e.mToVertIndex : mFromVertIndex < e.mFromVertIndex;
+		}
+
+		bool operator==(const MeshEdge& e) const
+		{
+			return mFromVertIndex == e.mFromVertIndex 
+				&& mToVertIndex == e.mToVertIndex;
+		}
+
+		int mFromVertIndex, mToVertIndex; 
+		int mHalfEdgeIndex; 
+	};
+
+	// check if the edge is following triangle order or not
+	bool checkEdgeOrientation(const MeshEdge &e, const shdfnd::Array<uint32_t> &indices)
+	{		
+		int offset0 = e.mHalfEdgeIndex % 3;
+		int offset1 = (offset0 < 2) ? 1 : -2;
+
+		int v0 = (int)indices[uint32_t(e.mHalfEdgeIndex)];
+		int v1 = (int)indices[uint32_t(e.mHalfEdgeIndex + offset1)];
+
+		if ((e.mFromVertIndex == v0) && (e.mToVertIndex == v1))
+			return true;
+
+		return false;
+	}
+
+	// check if two index pairs represent same edge regardless of order.
+	inline bool checkEdge(int ei0, int ei1, int ej0, int ej1)
+	{
+		return ( (ei0 == ej0) && (ei1 == ej1) ) ||
+			( (ei0 == ej1) && (ei1 == ej0) );
+	}
+
+	// compute ray edge intersection
+	bool intersectRayEdge(const PxVec3 &O, const PxVec3 &D, const PxVec3 &A, const PxVec3 &B, float &s, float &t)
+	{
+		// point on edge P = A + s * AB
+		// point on ray R = o + t * d
+		// for this two points to intersect, we have
+		// |AB -d| | s t | = o - A
+		const float eps = 1e-4;
+
+		PxVec3 OA = O - A;
+		PxVec3 AB = B - A;
+
+		float a = AB.dot(AB), b = -AB.dot(D);
+		float c = b, d = D.dot(D);
+
+		float e = AB.dot(OA);
+		float f = -D.dot(OA);
+
+		float det = a * d - b * c;
+		if (fabs(det) < eps) // coplanar case
+			return false;
+
+		float iPX_det = 1.0f / det;
+
+		s = (d * iPX_det) * e + (-b * iPX_det) * f;
+		t = (-c * iPX_det) * e + (a * iPX_det) * f;
+
+		return true;
+	}
+}
+
+
+struct nvidia::PxClothGeodesicTetherCookerImpl
+{
+
+	PxClothGeodesicTetherCookerImpl(const PxClothMeshDesc& desc);
+
+	uint32_t	getCookerStatus() const;
+	uint32_t	getNbTethersPerParticle() const;
+	void	getTetherData(uint32_t* userTetherAnchors, float* userTetherLengths) const;
+
+public:
+	// input
+	const PxClothMeshDesc&  mDesc;
+
+	// internal variables
+	uint32_t					mNumParticles;
+	shdfnd::Array<PxVec3>   mVertices;
+	shdfnd::Array<uint32_t>	mIndices;
+	shdfnd::Array<uint8_t>     mAttached;
+	shdfnd::Array<uint32_t>	mFirstVertTriAdj;
+	shdfnd::Array<uint32_t>	mVertTriAdjs;
+	shdfnd::Array<uint32_t>	mTriNeighbors; // needs changing for non-manifold support
+
+	// error status
+	uint32_t					mCookerStatus;
+
+	// output
+	shdfnd::Array<uint32_t>	mTetherAnchors;
+	shdfnd::Array<float>	mTetherLengths;
+
+protected:
+	void	createTetherData(const PxClothMeshDesc &desc);
+	int		computeVertexIntersection(uint32_t parent, uint32_t src, PathIntersection &path);
+	int		computeEdgeIntersection(uint32_t parent, uint32_t edge, float in_s, PathIntersection &path);
+	float	computeGeodesicDistance(uint32_t i, uint32_t parent, int &errorCode);
+	uint32_t	findTriNeighbors();
+	void	findVertTriNeighbors();
+
+private:
+	PxClothGeodesicTetherCookerImpl& operator=(const PxClothGeodesicTetherCookerImpl&);
+};
+
+PxClothGeodesicTetherCooker::PxClothGeodesicTetherCooker(const PxClothMeshDesc& desc)
+: mImpl(new PxClothGeodesicTetherCookerImpl(desc))
+{
+}
+
+PxClothGeodesicTetherCooker::~PxClothGeodesicTetherCooker()
+{
+	delete mImpl;
+}
+
+uint32_t PxClothGeodesicTetherCooker::getCookerStatus() const
+{
+	return mImpl->getCookerStatus();
+}
+
+uint32_t PxClothGeodesicTetherCooker::getNbTethersPerParticle() const
+{
+	return mImpl->getNbTethersPerParticle();
+}
+
+void PxClothGeodesicTetherCooker::getTetherData(uint32_t* userTetherAnchors, float* userTetherLengths) const
+{
+	mImpl->getTetherData(userTetherAnchors, userTetherLengths);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+PxClothGeodesicTetherCookerImpl::PxClothGeodesicTetherCookerImpl(const PxClothMeshDesc &desc)
+	:mDesc(desc),
+	mCookerStatus(0)
+{
+	createTetherData(desc);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+void PxClothGeodesicTetherCookerImpl::createTetherData(const PxClothMeshDesc &desc)
+{
+	mNumParticles = desc.points.count;
+	
+	if (!desc.invMasses.data)
+		return;
+
+	// assemble points
+	mVertices.resize(mNumParticles);
+	mAttached.resize(mNumParticles);
+	PxStrideIterator<const PxVec3> pIt((const PxVec3*)desc.points.data, desc.points.stride);
+	PxStrideIterator<const float> wIt((const float*)desc.invMasses.data, desc.invMasses.stride);
+	for(uint32_t i=0; i<mNumParticles; ++i)
+	{
+		mVertices[i] = *pIt++;
+		mAttached[i] = uint8_t(wIt.ptr() ? (*wIt++ == 0.0f) : 0);
+	}
+
+	// build triangle indices
+	if(desc.flags & PxMeshFlag::e16_BIT_INDICES)
+		gatherIndices<uint16_t>(mIndices, desc.triangles, desc.quads);
+	else
+		gatherIndices<uint32_t>(mIndices, desc.triangles, desc.quads);
+
+	// build vertex-triangle adjacencies
+	findVertTriNeighbors();
+
+	// build triangle-triangle adjacencies
+	mCookerStatus = findTriNeighbors();
+	if (mCookerStatus != 0)
+		return;
+
+	// build adjacent vertex list
+	shdfnd::Array<uint32_t> valency(mNumParticles+1, 0);
+	shdfnd::Array<uint32_t> adjacencies;
+	if(desc.flags & PxMeshFlag::e16_BIT_INDICES)
+		gatherAdjacencies<uint16_t>(valency, adjacencies, desc.triangles, desc.quads);
+	else
+		gatherAdjacencies<uint32_t>(valency, adjacencies, desc.triangles, desc.quads);
+
+	// build unique neighbors from adjacencies
+	shdfnd::Array<uint32_t> mark(valency.size(), 0);
+	shdfnd::Array<uint32_t> neighbors; neighbors.reserve(adjacencies.size());
+	for(uint32_t i=1, j=0; i<valency.size(); ++i)
+	{
+		for(; j<valency[i]; ++j)
+		{
+			uint32_t k = adjacencies[j];
+			if(mark[k] != i)
+			{
+				mark[k] = i;
+				neighbors.pushBack(k);
+			}
+		}
+		valency[i] = neighbors.size();
+	}
+
+	// create islands of attachment points
+	shdfnd::Array<uint32_t> vertexIsland(mNumParticles);
+	shdfnd::Array<VertexDistanceCount> vertexIslandHeap;
+
+	// put all the attachments in heap
+	for (uint32_t i = 0; i < mNumParticles; ++i)
+	{
+		// we put each attached point with large distance so that 
+		// we can prioritize things that are added during mesh traversal.
+		vertexIsland[i] = uint32_t(-1);
+		if (mAttached[i])
+			vertexIslandHeap.pushBack(VertexDistanceCount((int)i, FLT_MAX, 0));
+	}
+	uint32_t attachedCnt = vertexIslandHeap.size();
+
+	// no attached vertices
+	if (vertexIslandHeap.empty())
+		return;
+
+	// identify islands of attached vertices
+	shdfnd::Array<uint32_t> islandIndices;
+	shdfnd::Array<uint32_t> islandFirst;
+	uint32_t islandCnt = 0;
+	uint32_t islandIndexCnt = 0;
+
+	islandIndices.reserve(attachedCnt);
+	islandFirst.reserve(attachedCnt+1);
+
+	// while the island heap is not empty
+	while (!vertexIslandHeap.empty())
+	{
+		// pop vi from heap
+		VertexDistanceCount vi = popHeap(vertexIslandHeap);
+
+		// new cluster
+		if (vertexIsland[(uint32_t)vi.vertNr] == uint32_t(-1))
+		{
+			islandFirst.pushBack(islandIndexCnt++);
+			vertexIsland[(uint32_t)vi.vertNr] = islandCnt++;
+			vi.distance = 0;
+			islandIndices.pushBack((uint32_t)vi.vertNr);
+		}
+		
+		// for each adjacent vj that's not visited
+		const uint32_t begin = (uint32_t)valency[(uint32_t)vi.vertNr];
+		const uint32_t end = (uint32_t)valency[uint32_t(vi.vertNr + 1)];
+		for (uint32_t j = begin; j < end; ++j)
+		{
+			const uint32_t vj = neighbors[j];
+
+			// do not expand unattached vertices
+			if (!mAttached[vj])
+				continue; 
+
+			// already visited
+			if (vertexIsland[vj] != uint32_t(-1))
+				continue;
+
+			islandIndices.pushBack(vj);
+			islandIndexCnt++;
+			vertexIsland[vj] = vertexIsland[uint32_t(vi.vertNr)];
+			pushHeap(vertexIslandHeap, VertexDistanceCount((int)vj, vi.distance + 1.0f, 0));
+		}
+	}
+
+	islandFirst.pushBack(islandIndexCnt);
+
+	PX_ASSERT(islandCnt == (islandFirst.size() - 1));
+
+	/////////////////////////////////////////////////////////
+	uint32_t bufferSize = mNumParticles * islandCnt;
+	PX_ASSERT(bufferSize > 0);
+
+	shdfnd::Array<float> vertexDistanceBuffer(bufferSize, PX_MAX_F32);
+	shdfnd::Array<uint32_t> vertexParentBuffer(bufferSize, 0);
+	shdfnd::Array<VertexDistanceCount> vertexHeap;
+
+	// now process each island 
+	for (uint32_t i = 0; i < islandCnt; i++)
+	{
+		vertexHeap.clear();
+		float* vertexDistance = &vertexDistanceBuffer[0] + (i * mNumParticles);
+		uint32_t* vertexParent = &vertexParentBuffer[0] + (i * mNumParticles);
+
+		// initialize parent and distance
+		for (uint32_t j = 0; j < mNumParticles; ++j)
+		{
+			vertexParent[j] = j;
+			vertexDistance[j] = PX_MAX_F32;
+		}
+
+		// put all the attached vertices in this island to heap
+		const uint32_t beginIsland = islandFirst[i];
+		const uint32_t endIsland = islandFirst[i+1];
+		for (uint32_t j = beginIsland; j < endIsland; j++)
+		{
+			uint32_t vj = islandIndices[j];
+			vertexDistance[vj] = 0.0f;
+			vertexHeap.pushBack(VertexDistanceCount((int)vj, 0.0f, 0));
+		}
+
+		// no attached vertices in this island (error?)
+		PX_ASSERT(vertexHeap.empty() == false);
+		if (vertexHeap.empty())
+			continue;
+
+		// while heap is not empty
+		while (!vertexHeap.empty())
+		{
+			// pop vi from heap
+			VertexDistanceCount vi = popHeap(vertexHeap);
+
+			// obsolete entry ( we already found better distance)
+			if (vi.distance > vertexDistance[vi.vertNr])
+				continue;
+
+			// for each adjacent vj that's not visited
+			const int32_t begin = (int32_t)valency[(uint32_t)vi.vertNr];
+			const int32_t end = (int32_t)valency[uint32_t(vi.vertNr + 1)];
+			for (int32_t j = begin; j < end; ++j)
+			{
+				const int32_t vj = (int32_t)neighbors[(uint32_t)j];
+				PxVec3 edge = mVertices[(uint32_t)vj] - mVertices[(uint32_t)vi.vertNr];
+				const float edgeLength = edge.magnitude();
+				float newDistance = vi.distance + edgeLength;
+
+				if (newDistance < vertexDistance[vj])
+				{
+					vertexDistance[vj] = newDistance;
+					vertexParent[vj] = vertexParent[vi.vertNr];
+	
+					pushHeap(vertexHeap, VertexDistanceCount(vj, newDistance, 0));
+				}
+			}
+		}
+	}
+
+	const uint32_t maxTethersPerParticle = 4; // max tethers
+	const uint32_t nbTethersPerParticle = (islandCnt > maxTethersPerParticle) ? maxTethersPerParticle : islandCnt;
+
+	uint32_t nbTethers = nbTethersPerParticle * mNumParticles;
+	mTetherAnchors.resize(nbTethers);
+	mTetherLengths.resize(nbTethers);
+
+	// now process the parent and distance and add to fibers
+	for (uint32_t i = 0; i < mNumParticles; i++)
+	{
+		// we use the heap to sort out N-closest island
+		vertexHeap.clear();
+		for (uint32_t j = 0; j < islandCnt; j++)
+		{
+			int parent = (int)vertexParentBuffer[j * mNumParticles + i];
+			float edgeDistance = vertexDistanceBuffer[j * mNumParticles + i];
+			pushHeap(vertexHeap, VertexDistanceCount(parent, edgeDistance, 0));
+		}
+
+		// take out N-closest island from the heap
+		for (uint32_t j = 0; j < nbTethersPerParticle; j++)
+		{
+			VertexDistanceCount vi = popHeap(vertexHeap);
+			uint32_t parent = (uint32_t)vi.vertNr;
+			float distance = 0.0f;
+		
+			if (parent != i) 
+			{
+				float euclideanDistance = (mVertices[i] - mVertices[parent]).magnitude();
+				float dijkstraDistance = vi.distance;
+				int errorCode = 0;
+				float geodesicDistance = computeGeodesicDistance(i,parent, errorCode);
+				if (errorCode < 0)
+					geodesicDistance = dijkstraDistance;
+				distance = PxMax(euclideanDistance, geodesicDistance);
+			}
+
+			uint32_t tetherLoc = j * mNumParticles + i;
+			mTetherAnchors[ tetherLoc ] = parent;
+			mTetherLengths[ tetherLoc ] = distance;
+		}
+	}
+}
+
+///////////////////////////////////////////////////////////////////////////////
+uint32_t PxClothGeodesicTetherCookerImpl::getCookerStatus() const
+{
+	return mCookerStatus;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+uint32_t PxClothGeodesicTetherCookerImpl::getNbTethersPerParticle() const
+{
+	return mTetherAnchors.size() / mNumParticles;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+void  
+PxClothGeodesicTetherCookerImpl::getTetherData(uint32_t* userTetherAnchors, float* userTetherLengths) const
+{
+	intrinsics::memCopy(userTetherAnchors, mTetherAnchors.begin(), mTetherAnchors.size() * sizeof(uint32_t));
+	intrinsics::memCopy(userTetherLengths, mTetherLengths.begin(), mTetherLengths.size() * sizeof(float));
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// find triangle-triangle adjacency (return non-zero if there is an error)
+uint32_t PxClothGeodesicTetherCookerImpl::findTriNeighbors()
+{
+	shdfnd::Array<MeshEdge> edges;
+
+	mTriNeighbors.resize(mIndices.size(), uint32_t(-1));
+
+	// assemble all edges
+	uint32_t numTriangles = mIndices.size() / 3;
+	for (uint32_t i = 0; i < numTriangles; ++i)
+	{
+		uint32_t i0 = mIndices[3 * i];
+		uint32_t i1 = mIndices[3 * i + 1];
+		uint32_t i2 = mIndices[3 * i + 2];
+		edges.pushBack(MeshEdge((int)i0, (int)i1, int(3*i)));
+		edges.pushBack(MeshEdge((int)i1, (int)i2, int(3*i+1)));
+		edges.pushBack(MeshEdge((int)i2, (int)i0, int(3*i+2)));
+	}
+
+	shdfnd::sort(edges.begin(), edges.size());
+
+	int numEdges = (int)edges.size();
+	for(int i=0; i < numEdges; )
+	{
+		const MeshEdge& e0 = edges[(uint32_t)i];
+		bool orientation0 = checkEdgeOrientation(e0, mIndices);
+
+		int j = i;
+		while(++i < numEdges && edges[(uint32_t)i] == e0)
+			;
+
+		if(i - j > 2)
+			return 1; // non-manifold
+	
+		while(++j < i)
+		{
+			const MeshEdge& e1 = edges[(uint32_t)j];
+			bool orientation1 = checkEdgeOrientation(e1, mIndices);
+			mTriNeighbors[(uint32_t)e0.mHalfEdgeIndex] = (uint32_t)e1.mHalfEdgeIndex/3;
+			mTriNeighbors[(uint32_t)e1.mHalfEdgeIndex] = (uint32_t)e0.mHalfEdgeIndex/3;
+
+			if (orientation0 == orientation1)
+				return 2; // bad winding
+		}
+	}
+
+	return 0;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// find vertex triangle adjacency information
+void PxClothGeodesicTetherCookerImpl::findVertTriNeighbors()
+{
+	shdfnd::Array<VertTriangle> vertTriangles;
+	vertTriangles.reserve(mIndices.size());
+
+	int numTriangles = (int)mIndices.size() / 3;
+	for (int i = 0; i < numTriangles; ++i)
+	{
+		vertTriangles.pushBack(VertTriangle((int)mIndices[uint32_t(3*i)], i));
+		vertTriangles.pushBack(VertTriangle((int)mIndices[uint32_t(3*i+1)], i));
+		vertTriangles.pushBack(VertTriangle((int)mIndices[uint32_t(3*i+2)], i));
+	}
+
+	shdfnd::sort(vertTriangles.begin(), vertTriangles.size(), shdfnd::Less<VertTriangle>());
+	mFirstVertTriAdj.resize(mNumParticles);
+	mVertTriAdjs.reserve(mIndices.size());
+
+	for (uint32_t i = 0; i < (uint32_t)vertTriangles.size(); )
+	{
+		int v = vertTriangles[i].mVertIndex;
+
+		mFirstVertTriAdj[(uint32_t)v] = i;
+
+		while ((i < mIndices.size()) && (vertTriangles[i].mVertIndex == v))
+		{
+			int t = vertTriangles[i].mTriangleIndex;
+			mVertTriAdjs.pushBack((uint32_t)t);
+			i++;
+		}
+	}
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// compute intersection of a ray from a source vertex in direction toward parent
+int PxClothGeodesicTetherCookerImpl::computeVertexIntersection(uint32_t parent, uint32_t src, PathIntersection &path)
+{
+	if (src == parent)
+	{
+		path = PathIntersection(true, src, 0.0);
+		return 0;
+	}
+
+	float maxdot = -1.0f;
+	int closestVert = -1;
+
+	// gradient is toward the parent vertex
+	PxVec3 g = (mVertices[parent] - mVertices[src]).getNormalized();
+
+	// for every triangle incident on this vertex, we intersect against opposite edge of the triangle
+	uint32_t sfirst = mFirstVertTriAdj[src];
+	uint32_t slast = (src < ((uint32_t)mNumParticles-1)) ? mFirstVertTriAdj[src+1] : (uint32_t)mVertTriAdjs.size();
+	for (uint32_t adj = sfirst; adj < slast; adj++)
+	{
+		uint32_t tid = mVertTriAdjs[adj];
+		
+		uint32_t i0 = mIndices[tid*3];
+		uint32_t i1 = mIndices[tid*3+1];
+		uint32_t i2 = mIndices[tid*3+2];
+
+		int eid = 0;
+		if (i0 == src) eid = 1;
+		else if (i1 == src) eid = 2;
+		else if (i2 == src) eid = 0;
+		else continue; // error
+
+		// reshuffle so that src is located at i2
+		i0 = mIndices[tid*3 + eid];
+		i1 = mIndices[tid*3 + (eid+1)%3];
+		i2 = src;
+
+		PxVec3 p0 = mVertices[i0];
+		PxVec3 p1 = mVertices[i1];
+		PxVec3 p2 = mVertices[i2];
+
+		// check if we hit source immediately from this triangle
+		if (i0 == parent)
+		{
+			path = PathIntersection(true, parent, (p0 - p2).magnitude());
+			return 1;
+		}
+
+		if (i1 == parent)
+		{
+			path = PathIntersection(true, parent, (p1 - p2).magnitude());
+			return 1;
+		}
+
+		// ray direction is the gradient projected on the plane of this triangle
+		PxVec3 n = ((p0 - p2).cross(p1 - p2)).getNormalized();
+		PxVec3 d = (g - g.dot(n) * n).getNormalized();
+
+		// find intersection of ray (p2, d) against the edge (p0,p1)
+		float s, t;
+		bool result = intersectRayEdge(p2, d, p0, p1, s, t);
+		if (result == false)
+			continue;
+
+		// t should be positive, otherwise we just hit the triangle in opposite direction, so ignore
+		const float eps = 1e-5;
+		if (t > -eps)
+		{		
+			PxVec3 ip; // intersection point
+			if (( s > -eps ) && (s < (1.0f + eps)))
+			{				
+				// if intersection point is too close to each vertex, we record a vertex intersection
+				if ( ( s < eps) || (s > (1.0f-eps)))
+				{
+					path.vertOrTriangle = true;
+					path.index = (s < eps) ? i0 : i1;
+					path.distance = (p2 - mVertices[path.index]).magnitude();				
+				}
+				else // found an edge instersection
+				{
+					ip = p0 + s * (p1 - p0);
+					path = PathIntersection(false, tid*3 + eid, (p2 - ip).magnitude(), s);					
+				}
+				return 1;
+			}
+		}
+
+		// for fall back (see below)
+		PxVec3 d0 = (p0 - p2).getNormalized();
+		PxVec3 d1 = (p1 - p2).getNormalized();
+		float d0dotg = d0.dot(d);
+		float d1dotg = d1.dot(d);
+
+		if (d0dotg > maxdot)
+		{
+			closestVert = (int)i0;
+			maxdot = d0dotg;
+		}
+		if (d1dotg > maxdot)
+		{
+			closestVert = (int)i1;
+			maxdot = d1dotg;
+		}
+	} // end for (uint32_t adj = sfirst...
+
+	// Fall back to use greedy (Dijkstra-like) path selection. 
+	// This happens as triangles are curved and we may not find intersection on any triangle.
+	// In this case, we choose a vertex closest to the gradient direction.
+	if (closestVert > 0)
+	{
+		path = PathIntersection(true, (uint32_t)closestVert, (mVertices[src] - mVertices[(uint32_t)closestVert]).magnitude());
+		return 1;
+	}
+
+	// Error, (possibly dangling vertex)
+	return -1;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// compute intersection of a ray from a source vertex in direction toward parent
+int PxClothGeodesicTetherCookerImpl::computeEdgeIntersection(uint32_t parent, uint32_t edge, float in_s, PathIntersection &path)
+{
+	int tid = (int)edge / 3;
+	int eid = (int)edge % 3;
+
+	uint32_t e0 = mIndices[uint32_t(tid*3 + eid)];
+	uint32_t e1 = mIndices[uint32_t(tid*3 + (eid+1)%3)];
+
+	PxVec3 v0 = mVertices[e0];
+	PxVec3 v1 = mVertices[e1];
+
+	PxVec3 v = v0 + in_s * (v1 - v0);
+	PxVec3 g = mVertices[parent] - v;
+
+	uint32_t triNbr = mTriNeighbors[edge];
+
+	if (triNbr == uint32_t(-1)) // boundary edge
+	{
+		float dir = g.dot(v1-v0);
+		uint32_t vid = (dir > 0) ? e1 : e0;
+		path = PathIntersection(true, vid, (mVertices[vid] - v).magnitude());
+		return 1;
+	}
+
+	uint32_t i0 = mIndices[triNbr*3];
+	uint32_t i1 = mIndices[triNbr*3+1];
+	uint32_t i2 = mIndices[triNbr*3+2];
+
+	// vertex is sorted s.t i0,i1 contains the edge point
+	if ( checkEdge((int)i0, (int)i1, (int)e0, (int)e1)) {
+		eid = 0;
+	}
+	else if ( checkEdge((int)i1, (int)i2, (int)e0, (int)e1)) {
+		eid = 1;
+		uint32_t tmp = i2;
+		i2 = i0;
+		i0 = i1;
+		i1 = tmp;
+	}
+	else if ( checkEdge((int)i2, (int)i0, (int)e0, (int)e1)) 
+	{
+		eid = 2;
+		uint32_t tmp = i0;
+		i0 = i2;
+		i2 = i1;
+		i1 = tmp;
+	}
+
+	// we hit the parent
+	if (i2 == parent)
+	{
+		path = PathIntersection(true, i2, (mVertices[i2] - v).magnitude());
+		return 1;
+	}
+
+	PxVec3 p0 = mVertices[i0];
+	PxVec3 p1 = mVertices[i1];
+	PxVec3 p2 = mVertices[i2];
+
+	// project gradient vector on the plane of the triangle
+	PxVec3 n = ((p0 - p2).cross(p1 - p2)).getNormalized();
+	g = (g - g.dot(n) * n).getNormalized();
+
+	float s = 0.0f, t = 0.0f;
+	const float eps = 1e-5;
+	PxVec3 ip;
+
+	// intersect against edge form p2 to p0
+	if (intersectRayEdge(v, g, p2, p0, s, t) && ( s >= -eps) && ( s <= (1.0f+eps) ) && (t > -eps))
+	{
+		if ( ( s < eps) || (s > (1.0f-eps)))
+		{
+			path.vertOrTriangle = true;
+			path.index = (s < eps) ? i2 : i0;
+			path.distance = (mVertices[path.index] - v).magnitude();
+		}
+		else
+		{
+			ip = p2 + s * (p0 - p2);
+			path = PathIntersection(false, triNbr*3 + (eid + 2) % 3, (ip - v).magnitude(), s);
+			
+		}
+
+		return 1;
+	}
+
+	// intersect against edge form p1 to p2
+	if (intersectRayEdge(v, g, p1, p2, s, t) && ( s >= -eps) && ( s <= (1.0f+eps) ) && (t > -eps))
+	{
+		if ( ( s < eps) || (s > (1.0f-eps)))
+		{
+			path.vertOrTriangle = true;
+			path.index = (s < eps) ? i1 : i2;
+			path.distance = (mVertices[path.index] - v).magnitude();
+		}
+		else
+		{
+			ip = p1 + s * (p2 - p1);
+			path = PathIntersection(false, triNbr*3 + (eid + 1) % 3, (ip - v).magnitude(), s);
+		}
+		
+		return 1;
+	}
+
+	// fallback to pick closer vertex when no edges intersect
+	float dir = g.dot(v1-v0);
+	path.vertOrTriangle = true;
+	path.index = (dir > 0) ? e1 : e0;
+	path.distance = (mVertices[path.index] - v).magnitude();
+
+	return 1;
+}
+
+
+///////////////////////////////////////////////////////////////////////////////
+// compute geodesic distance and path from vertex i to its parent
+float PxClothGeodesicTetherCookerImpl::computeGeodesicDistance(uint32_t i, uint32_t parent, int &errorCode)
+{
+	if (i == parent)
+		return 0.0f;
+		
+	PathIntersection path;
+	
+	errorCode = 0;
+
+	// find intial intersection
+	int status = computeVertexIntersection(parent, i, path);
+	if (status < 0)
+	{
+		errorCode = -1;
+		return 0;
+	}
+
+	int pathcnt = 0;
+	float geodesicDistance = 0;
+
+	while (status > 0)
+	{	
+		geodesicDistance += path.distance;
+
+		if (path.vertOrTriangle)
+			status = computeVertexIntersection(parent, path.index, path);
+		else 
+			status = computeEdgeIntersection(parent, path.index, path.s, path);
+
+		// cannot find valid path
+		if (status < 0) 
+		{
+			errorCode = -2;
+			return 0.0f;
+		}
+
+		// possibly cycles, too many path
+		if (pathcnt > 1000) 
+		{		
+			errorCode = -3;
+			return 0.0f;
+		}
+
+		pathcnt++;
+	}
+
+	return geodesicDistance;
+}
+
+
+
+
+#endif //APEX_USE_CLOTH_API
+
+