diff options
Diffstat (limited to 'Runtime/Geometry')
24 files changed, 6097 insertions, 0 deletions
diff --git a/Runtime/Geometry/AABB.cpp b/Runtime/Geometry/AABB.cpp new file mode 100644 index 0000000..c8e01b4 --- /dev/null +++ b/Runtime/Geometry/AABB.cpp @@ -0,0 +1,241 @@ +#include "UnityPrefix.h" +#include "AABB.h" +#include "Runtime/Math/Quaternion.h" + +const AABB AABB::zero = AABB(Vector3f::zero, Vector3f::zero); + +void CalculateClosestPoint (const Vector3f& rkPoint, const AABB& rkBox, Vector3f& outPoint, float& outSqrDistance) +{ + // compute coordinates of point in box coordinate system + Vector3f kClosest = rkPoint - rkBox.GetCenter(); + + // project test point onto box + float fSqrDistance = 0.0f; + float fDelta; + + for (int i=0;i<3;i++) + { + if ( kClosest[i] < -rkBox.GetExtent (i) ) + { + fDelta = kClosest[i] + rkBox.GetExtent (i); + fSqrDistance += fDelta * fDelta; + kClosest[i] = -rkBox.GetExtent (i); + } + else if ( kClosest[i] > rkBox.GetExtent(i) ) + { + fDelta = kClosest[i] - rkBox.GetExtent (i); + fSqrDistance += fDelta * fDelta; + kClosest[i] = rkBox.GetExtent (i); + } + } + + // Inside + if (fSqrDistance == 0.0F) + { + outPoint = rkPoint; + outSqrDistance = 0.0F; + } + // Outside + else + { + outPoint = kClosest + rkBox.GetCenter(); + outSqrDistance = fSqrDistance; + } +} + + +// Sphere-AABB distance, Arvo's algorithm +float CalculateSqrDistance (const Vector3f& rkPoint, const AABB& rkBox) +{ + Vector3f closest = rkPoint - rkBox.GetCenter(); + float sqrDistance = 0.0f; + + for (int i = 0; i < 3; ++i) + { + float clos = closest[i]; + float ext = rkBox.GetExtent(i); + if (clos < -ext) + { + float delta = clos + ext; + sqrDistance += delta * delta; + closest[i] = -ext; + } + else if (clos > ext) + { + float delta = clos - ext; + sqrDistance += delta * delta; + closest[i] = ext; + } + } + + return sqrDistance; +} + +void AABB::GetVertices (Vector3f* outVertices) const +{ + outVertices[0] = m_Center + Vector3f (-m_Extent.x, -m_Extent.y, -m_Extent.z); + outVertices[1] = m_Center + Vector3f (+m_Extent.x, -m_Extent.y, -m_Extent.z); + outVertices[2] = m_Center + Vector3f (-m_Extent.x, +m_Extent.y, -m_Extent.z); + outVertices[3] = m_Center + Vector3f (+m_Extent.x, +m_Extent.y, -m_Extent.z); + + outVertices[4] = m_Center + Vector3f (-m_Extent.x, -m_Extent.y, +m_Extent.z); + outVertices[5] = m_Center + Vector3f (+m_Extent.x, -m_Extent.y, +m_Extent.z); + outVertices[6] = m_Center + Vector3f (-m_Extent.x, +m_Extent.y, +m_Extent.z); + outVertices[7] = m_Center + Vector3f (+m_Extent.x, +m_Extent.y, +m_Extent.z); +} + +void MinMaxAABB::GetVertices( Vector3f outVertices[8] ) const +{ + // 7-----6 + // / /| + // 3-----2 | + // | 4 | 5 + // | |/ + // 0-----1 + outVertices[0].Set( m_Min.x, m_Min.y, m_Min.z ); + outVertices[1].Set( m_Max.x, m_Min.y, m_Min.z ); + outVertices[2].Set( m_Max.x, m_Max.y, m_Min.z ); + outVertices[3].Set( m_Min.x, m_Max.y, m_Min.z ); + outVertices[4].Set( m_Min.x, m_Min.y, m_Max.z ); + outVertices[5].Set( m_Max.x, m_Min.y, m_Max.z ); + outVertices[6].Set( m_Max.x, m_Max.y, m_Max.z ); + outVertices[7].Set( m_Min.x, m_Max.y, m_Max.z ); +} + + + +bool AABB::IsInside (const Vector3f& inPoint) const +{ + if (inPoint[0] < m_Center[0] - m_Extent[0]) + return false; + if (inPoint[0] > m_Center[0] + m_Extent[0]) + return false; + + if (inPoint[1] < m_Center[1] - m_Extent[1]) + return false; + if (inPoint[1] > m_Center[1] + m_Extent[1]) + return false; + + if (inPoint[2] < m_Center[2] - m_Extent[2]) + return false; + if (inPoint[2] > m_Center[2] + m_Extent[2]) + return false; + + return true; +} + +void AABB::Encapsulate (const Vector3f& inPoint) { + MinMaxAABB temp = *this; + temp.Encapsulate (inPoint); + FromMinMaxAABB (temp); +} + +bool MinMaxAABB::IsInside (const Vector3f& inPoint) const +{ + if (inPoint[0] < m_Min[0]) + return false; + if (inPoint[0] > m_Max[0]) + return false; + + if (inPoint[1] < m_Min[1]) + return false; + if (inPoint[1] > m_Max[1]) + return false; + + if (inPoint[2] < m_Min[2]) + return false; + if (inPoint[2] > m_Max[2]) + return false; + + return true; +} + +MinMaxAABB AddAABB (const MinMaxAABB& lhs, const MinMaxAABB& rhs) +{ + MinMaxAABB minMax; + if (lhs.IsValid()) + minMax = lhs; + + if (rhs.IsValid()) + { + minMax.Encapsulate (rhs.GetMax ()); + minMax.Encapsulate (rhs.GetMin ()); + } + + return minMax; +} + +inline Vector3f RotateExtents (const Vector3f& extents, const Matrix3x3f& rotation) +{ + Vector3f newExtents; + for (int i=0;i<3;i++) + newExtents[i] = Abs (rotation.Get (i, 0) * extents.x) + Abs (rotation.Get (i, 1) * extents.y) + Abs (rotation.Get (i, 2) * extents.z); + return newExtents; +} + +inline Vector3f RotateExtents (const Vector3f& extents, const Matrix4x4f& rotation) +{ + Vector3f newExtents; + for (int i=0;i<3;i++) + newExtents[i] = Abs (rotation.Get (i, 0) * extents.x) + Abs (rotation.Get (i, 1) * extents.y) + Abs (rotation.Get (i, 2) * extents.z); + return newExtents; +} + +void TransformAABB (const AABB& aabb, const Vector3f& position, const Quaternionf& rotation, AABB& result) +{ + Matrix3x3f m; + QuaternionToMatrix (rotation, m); + + Vector3f extents = RotateExtents (aabb.GetExtent (), m); + Vector3f center = m.MultiplyPoint3 (aabb.GetCenter ()); + center += position; + result.SetCenterAndExtent( center, extents ); +} + +void TransformAABB (const AABB& aabb, const Matrix4x4f& transform, AABB& result) +{ + Vector3f extents = RotateExtents (aabb.GetExtent (), transform); + Vector3f center = transform.MultiplyPoint3 (aabb.GetCenter ()); + result.SetCenterAndExtent( center, extents ); +} + +void TransformAABBSlow (const AABB& aabb, const Matrix4x4f& transform, AABB& result) +{ + MinMaxAABB transformed; + transformed.Init (); + + Vector3f v[8]; + aabb.GetVertices (v); + for (int i=0;i<8;i++) + transformed.Encapsulate (transform.MultiplyPoint3 (v[i])); + + result = transformed; +} + + +void InverseTransformAABB (const AABB& aabb, const Vector3f& position, const Quaternionf& rotation, AABB& result) +{ + Matrix3x3f m; + QuaternionToMatrix (Inverse (rotation), m); + + Vector3f extents = RotateExtents (aabb.GetExtent (), m); + Vector3f center = aabb.GetCenter () - position; + center = m.MultiplyPoint3 (center); + + result.SetCenterAndExtent( center, extents ); +} + +bool IsContainedInAABB (const AABB& inside, const AABB& bigBounds) +{ + bool outside = false; + outside |= inside.m_Center[0] - inside.m_Extent[0] < bigBounds.m_Center[0] - bigBounds.m_Extent[0]; + outside |= inside.m_Center[0] + inside.m_Extent[0] > bigBounds.m_Center[0] + bigBounds.m_Extent[0]; + + outside |= inside.m_Center[1] - inside.m_Extent[1] < bigBounds.m_Center[1] - bigBounds.m_Extent[1]; + outside |= inside.m_Center[1] + inside.m_Extent[1] > bigBounds.m_Center[1] + bigBounds.m_Extent[1]; + + outside |= inside.m_Center[2] - inside.m_Extent[2] < bigBounds.m_Center[2] - bigBounds.m_Extent[2]; + outside |= inside.m_Center[2] + inside.m_Extent[2] > bigBounds.m_Center[2] + bigBounds.m_Extent[2]; + + return !outside; +} diff --git a/Runtime/Geometry/AABB.h b/Runtime/Geometry/AABB.h new file mode 100644 index 0000000..8b80cf1 --- /dev/null +++ b/Runtime/Geometry/AABB.h @@ -0,0 +1,202 @@ +#ifndef AABB_H +#define AABB_H + +#include "Runtime/Math/Vector3.h" +#include "Runtime/Math/Matrix4x4.h" + +class MinMaxAABB; +class Quaternionf; +class Matrix3x3f; + + +class AABB +{ +public: + static const AABB zero; + +public: + Vector3f m_Center; + Vector3f m_Extent; + + + AABB () {} + AABB (const Vector3f& c, const Vector3f& e) { m_Center = c; m_Extent = e; } + DECLARE_SERIALIZE_OPTIMIZE_TRANSFER (AABB) + + AABB (const MinMaxAABB& aabb) { FromMinMaxAABB (aabb); } + AABB& operator = (const MinMaxAABB& aabb) { FromMinMaxAABB (aabb); return *this; } + + bool operator == (const AABB& b) const { return m_Center == b.m_Center && m_Extent == b.m_Extent; } + + void SetCenterAndExtent( const Vector3f& c, const Vector3f& e ) { m_Center = c; m_Extent = e; } + + Vector3f& GetCenter () { return m_Center; } + Vector3f& GetExtent () { return m_Extent; } + float& GetExtent (int i) { return m_Extent[i]; } + const Vector3f& GetCenter ()const { return m_Center; } + const Vector3f& GetExtent ()const { return m_Extent; } + float GetExtent (int i)const { return m_Extent[i]; } + + Vector3f GetMin () const { return m_Center - m_Extent; } + Vector3f GetMax () const { return m_Center + m_Extent; } + + void Expand (float inValue); + + bool IsValid () const; + + bool IsInside (const Vector3f& inPoint) const; + + void GetVertices (Vector3f* outVertices) const; + void EXPORT_COREMODULE Encapsulate (const Vector3f& inPoint); + + private: + + void FromMinMaxAABB (const MinMaxAABB& aabb); +}; + +bool IsContainedInAABB (const AABB& inside, const AABB& bigBounds); + + +// Find minimum AABB which includes both AABBs +MinMaxAABB AddAABB (const MinMaxAABB& lhs, const MinMaxAABB& rhs); + +// Transforms AABB. +// Can be thought of as Converting OBB to an AABB: +// rotate the center and extents of the OBB And find the smallest enclosing AABB around it. +void TransformAABB (const AABB& aabb, const Vector3f& position, const Quaternionf& rotation, AABB& result); + +/// This is not mathematically correct for non-uniform scaled objects. But it seems to work well enough. +/// If you use it with non-uniform scale make sure to verify it extensively. +EXPORT_COREMODULE void TransformAABB (const AABB& aabb, const Matrix4x4f& transform, AABB& result); + +/// This version is much slower but works correctly with non-uniform scale +void TransformAABBSlow (const AABB& aabb, const Matrix4x4f& transform, AABB& result); + +void InverseTransformAABB (const AABB& aabb, const Vector3f& position, const Quaternionf& rotation, AABB& result); + + +/// The closest distance to the surface or inside the aabb. +float CalculateSqrDistance (const Vector3f& rkPoint, const AABB& rkBox); + + +/// Returns the sqr distance and the closest point inside or on the surface of the aabb. +/// If inside the aabb, distance will be zero and rkPoint will be returned. +EXPORT_COREMODULE void CalculateClosestPoint (const Vector3f& rkPoint, const AABB& rkBox, Vector3f& outPoint, float& outSqrDistance); + +class MinMaxAABB +{ +public: + + Vector3f m_Min; + Vector3f m_Max; + + MinMaxAABB () { Init (); } + MinMaxAABB (Vector3f min, Vector3f max) : m_Min(min), m_Max(max) { }; + MinMaxAABB (const AABB& aabb) { FromAABB (aabb); } + MinMaxAABB& operator = (const AABB& aabb) { FromAABB (aabb); return *this; } + //DECLARE_SERIALIZE_OPTIMIZE_TRANSFER (MinMaxAABB) + + void Init (); + + const Vector3f& GetMin () const { return m_Min; } + const Vector3f& GetMax () const { return m_Max; } + Vector3f GetCenter () const { return 0.5F * (m_Max + m_Min); } + Vector3f GetExtent () const { return 0.5F * (m_Max - m_Min); } + Vector3f GetSize () const { return (m_Max - m_Min); } + + void Encapsulate (const Vector3f& inPoint); + void Encapsulate (const AABB& aabb); + void Encapsulate (const MinMaxAABB& other); + + void Expand (float inValue); + void Expand (const Vector3f& inOffset); + + // TODO : rename - it has different meaning than AABB::IsValid + bool IsValid () const; + + bool IsInside (const Vector3f& inPoint) const; + + void GetVertices( Vector3f outVertices[8] ) const; + +private: + + void FromAABB (const AABB& inAABB); +}; + +inline void AABB::Expand (float inValue) +{ + m_Extent += Vector3f (inValue, inValue, inValue); +} + +inline void AABB::FromMinMaxAABB (const MinMaxAABB& inAABB) +{ + m_Center = (inAABB.GetMax () + inAABB.GetMin ()) * 0.5F; + m_Extent = (inAABB.GetMax () - inAABB.GetMin ()) * 0.5F; +} + +inline bool AABB::IsValid () const +{ + return IsFinite(m_Center) && IsFinite(m_Extent); +} + +inline void MinMaxAABB::Encapsulate (const Vector3f& inPoint) +{ + m_Min = min (m_Min, inPoint); + m_Max = max (m_Max, inPoint); +} + +inline void MinMaxAABB::Encapsulate (const AABB& aabb) +{ + Encapsulate (aabb.GetCenter()+aabb.GetExtent()); + Encapsulate (aabb.GetCenter()-aabb.GetExtent()); +} + +inline void MinMaxAABB::Encapsulate (const MinMaxAABB& other) +{ + m_Min = min (m_Min, other.m_Min); + m_Max = max (m_Max, other.m_Max); +} + +inline void MinMaxAABB::Expand (float inValue) +{ + Vector3f offset = Vector3f (inValue, inValue, inValue); + m_Min -= offset; + m_Max += offset; +} + +inline void MinMaxAABB::Expand (const Vector3f& inOffset) +{ + m_Min -= inOffset; + m_Max += inOffset; +} + +inline bool MinMaxAABB::IsValid () const +{ + return !(m_Min == Vector3f::infinityVec || m_Max == -Vector3f::infinityVec); +} + +inline void MinMaxAABB::Init () +{ + m_Min = Vector3f::infinityVec; + m_Max = -Vector3f::infinityVec; +} + +inline void MinMaxAABB::FromAABB (const AABB& inAABB) +{ + m_Min = inAABB.GetCenter () - inAABB.GetExtent (); + m_Max = inAABB.GetCenter () + inAABB.GetExtent (); +} +template<class TransferFunction> inline +void AABB::Transfer (TransferFunction& transfer) +{ + TRANSFER (m_Center); + TRANSFER (m_Extent); +} +/*template<class TransferFunction> inline +void MinMaxAABB::Transfer (TransferFunction& transfer) +{ + TRANSFER (m_Min); + TRANSFER (m_Max); +}*/ + +#endif diff --git a/Runtime/Geometry/BoundingUtils.cpp b/Runtime/Geometry/BoundingUtils.cpp new file mode 100644 index 0000000..7c9b361 --- /dev/null +++ b/Runtime/Geometry/BoundingUtils.cpp @@ -0,0 +1,443 @@ +#include "UnityPrefix.h" +#include "BoundingUtils.h" +#include "Plane.h" +#include "AABB.h" +#include "Intersection.h" + +// -------------------------------------------------------------------------- + +void GetFrustumPoints( const Matrix4x4f& clipToWorld, Vector3f* frustum ) +{ + clipToWorld.PerspectiveMultiplyPoint3( Vector3f(-1,-1,-1), frustum[0] ); + clipToWorld.PerspectiveMultiplyPoint3( Vector3f( 1,-1,-1), frustum[1] ); + clipToWorld.PerspectiveMultiplyPoint3( Vector3f( 1, 1,-1), frustum[2] ); + clipToWorld.PerspectiveMultiplyPoint3( Vector3f(-1, 1,-1), frustum[3] ); + clipToWorld.PerspectiveMultiplyPoint3( Vector3f(-1,-1, 1), frustum[4] ); + clipToWorld.PerspectiveMultiplyPoint3( Vector3f( 1,-1, 1), frustum[5] ); + clipToWorld.PerspectiveMultiplyPoint3( Vector3f( 1, 1, 1), frustum[6] ); + clipToWorld.PerspectiveMultiplyPoint3( Vector3f(-1, 1, 1), frustum[7] ); +} + +void GetFrustumPortion( const Vector3f* frustum, float nearSplit, float farSplit, Vector3f* outPortion ) +{ + outPortion[0] = Lerp( frustum[0], frustum[0+4], nearSplit ); + outPortion[1] = Lerp( frustum[1], frustum[1+4], nearSplit ); + outPortion[2] = Lerp( frustum[2], frustum[2+4], nearSplit ); + outPortion[3] = Lerp( frustum[3], frustum[3+4], nearSplit ); + outPortion[4] = Lerp( frustum[0], frustum[0+4], farSplit ); + outPortion[5] = Lerp( frustum[1], frustum[1+4], farSplit ); + outPortion[6] = Lerp( frustum[2], frustum[2+4], farSplit ); + outPortion[7] = Lerp( frustum[3], frustum[3+4], farSplit ); +} + +static bool ClipTest( const float p, const float q, float& u1, float& u2 ) +{ + // Return value is 'true' if line segment intersects the current test + // plane. Otherwise 'false' is returned in which case the line segment + // is entirely clipped. + const float EPS = 1.0e-10f; + if( p < -EPS ) { + float r = q/p; + if( r > u2 ) + return false; + else { + if( r > u1 ) + u1 = r; + return true; + } + } + else if( p > EPS ) + { + float r = q/p; + if( r < u1 ) + return false; + else { + if( r < u2 ) + u2 = r; + return true; + } + } + else + { + return q >= 0.0f; + } +} + +static bool IntersectLineAABB( Vector3f& v, const Vector3f& p, const Vector3f& dir, const MinMaxAABB& b ) +{ + float t1 = 0.0f; + float t2 = 1.0e30f; + bool intersect = false; + if (ClipTest(-dir.z,p.z-b.GetMin().z,t1,t2) && ClipTest(dir.z,b.GetMax().z-p.z,t1,t2) && + ClipTest(-dir.y,p.y-b.GetMin().y,t1,t2) && ClipTest(dir.y,b.GetMax().y-p.y,t1,t2) && + ClipTest(-dir.x,p.x-b.GetMin().x,t1,t2) && ClipTest(dir.x,b.GetMax().x-p.x,t1,t2)) + { + if( 0 <= t1 ) { + v = p + dir * t1; + intersect = true; + } + if( 0 <= t2 ) { + v = p + dir * t2; + intersect = true; + } + } + return intersect; +} + +inline bool ClipPolysByPlane( UInt8 numPoints, const Vector3f* __restrict input, const Plane& A, UInt8* __restrict pNumOutPoints, Vector3f* __restrict output, UInt8* __restrict pNumIntermPoints, Vector3f* __restrict interm ) +{ + int i; + if( numPoints < 3 ) + { + *pNumOutPoints = 0; + return false; + } + + UInt8 numOutPoints = 0; + UInt8& numIntermPoints = *pNumIntermPoints; + + Vector3f temp; + + bool * outside = (bool*)alloca(numPoints*sizeof(bool)); + for( i = 0; i < numPoints; ++i ) + outside[i] = A.GetDistanceToPoint( input[i] ) < 0.0f; + + for( i = 0; i < numPoints; ++i ) + { + int idNext = (i+1) % numPoints; + + // both outside -> save none + if(outside[i] && outside[idNext]) + continue; + + // outside-inside -> save intersection and i+1 + if(outside[i]) + { + if( IntersectSegmentPlane( input[i], input[idNext], A, &temp ) ) + { + output[numOutPoints++] = temp; + interm[numIntermPoints++] = temp; + } + + output[numOutPoints++] = input[idNext]; + continue; + } + + // inside-outside -> save intersection + if(outside[idNext]) + { + if( IntersectSegmentPlane( input[i], input[idNext], A, &temp ) ) + { + output[numOutPoints++] = temp; + interm[numIntermPoints++] = temp; + } + + continue; + } + + output[numOutPoints++] = input[idNext]; + } + *pNumOutPoints = numOutPoints; + return numOutPoints ? true : false; +} + +void CalcHullBounds(const Vector3f* __restrict hullPoints, const UInt8* __restrict hullCounts, UInt8 hullFaces, const Plane& nearPlane, const Matrix4x4f& cameraWorldToClip, MinMaxAABB& aabb) +{ + Vector3f outputData[128]; + Vector3f intermData[128]; + + UInt8 outputDataCounts[64]; + UInt8 intermDataCounts[64]; + + const Vector3f* __restrict inputPoints = hullPoints; + const UInt8* __restrict inputCounts = hullCounts; + + Vector3f* __restrict outputPoints = outputData; + UInt8* __restrict outputCounts = outputDataCounts; + intermDataCounts[0] = 0; + for( UInt8 i = 0; i != hullFaces; ++i ) + { + const UInt8 inputCount = *inputCounts++; + ClipPolysByPlane( inputCount, inputPoints, nearPlane, outputCounts, outputPoints, intermDataCounts, intermData); + inputPoints += inputCount; + outputPoints += *outputCounts; + outputCounts++; + } + + // Project hull's points onto the screen and compute bounding rectangle of them. + Vector3f projectedPoint; + outputPoints = outputData; + outputCounts = outputDataCounts; + + for( int f = 0; f < hullFaces; ++f ) + { + const UInt8 numPoints = *outputCounts++; + for( int i = 0; i < numPoints; ++i ) + { + cameraWorldToClip.PerspectiveMultiplyPoint3( outputPoints[i], projectedPoint ); + aabb.Encapsulate( projectedPoint ); + } + outputPoints += numPoints; + } + + if( aabb.m_Min.x < -1.0f ) aabb.m_Min.x = -1.0f; + if( aabb.m_Min.y < -1.0f ) aabb.m_Min.y = -1.0f; + if( aabb.m_Max.x > 1.0f ) aabb.m_Max.x = 1.0f; + if( aabb.m_Max.y > 1.0f ) aabb.m_Max.y = 1.0f; + +} + + +// Returns our "focused region of interest": frustum, clipped by scene bounds, +// extruded towards light and clipped by scene bounds again. +// +// Frustum is +// { -1, -1, -1 } +// { 1, -1, -1 } +// { 1, 1, -1 } +// { -1, 1, -1 } +// { -1, -1, 1 } +// { 1, -1, 1 } +// { 1, 1, 1 } +// { -1, 1, 1 } + + + +void CalculateFocusedLightHull( const Vector3f* frustum, const Vector3f& lightDir, const MinMaxAABB& sceneAABB, std::vector<Vector3f>& points ) +{ + int i; + Vector3f tempPoints[3][256]; + UInt8 tempCounts[3][128]; + + Plane planes[6]; + planes[0].SetABCD( 0, 1, 0, -sceneAABB.GetMin().y ); + planes[1].SetABCD( 0, -1, 0, sceneAABB.GetMax().y ); + planes[2].SetABCD( 1, 0, 0, -sceneAABB.GetMin().x ); + planes[3].SetABCD( -1, 0, 0, sceneAABB.GetMax().x ); + planes[4].SetABCD( 0, 0, 1, -sceneAABB.GetMin().z ); + planes[5].SetABCD( 0, 0, -1, sceneAABB.GetMax().z ); + + Vector3f* __restrict pData[2] = { (Vector3f*)&tempPoints[0][0], (Vector3f*)&tempPoints[1][0] }; + UInt8* __restrict pDataCounts[2] = { (UInt8*)&tempCounts[0][0], (UInt8*)&tempCounts[1][0] }; + + Vector3f* v = *pData; + UInt8* c = *pDataCounts; + + UInt32 numFaces = 6; + + c[0] = c[1] = c[2] = c[3] = c[4] = c[5] = 4; + v[ 0] = frustum[0]; v[ 1] = frustum[1]; v[ 2] = frustum[2]; v[ 3] = frustum[3]; + v[ 4] = frustum[7]; v[ 5] = frustum[6]; v[ 6] = frustum[5]; v[ 7] = frustum[4]; + v[ 8] = frustum[0]; v[ 9] = frustum[3]; v[10] = frustum[7]; v[11] = frustum[4]; + v[12] = frustum[1]; v[13] = frustum[5]; v[14] = frustum[6]; v[15] = frustum[2]; + v[16] = frustum[4]; v[17] = frustum[5]; v[18] = frustum[1]; v[19] = frustum[0]; + v[20] = frustum[6]; v[21] = frustum[7]; v[22] = frustum[3]; v[23] = frustum[2]; + + + UInt32 numTotalPoints; + UInt32 vIn = 0; + + + Vector3f* intermPoints = &tempPoints[2][0]; + UInt8* intermCounts = &tempCounts[2][0]; + + for( int p = 0; p < 6; ++p ) + { + const Vector3f* __restrict inputPoints=pData[vIn]; + Vector3f* __restrict outputPoints = pData[1-vIn]; + + UInt8* __restrict inputCounts=pDataCounts[vIn]; + UInt8* __restrict outputCounts = pDataCounts[1-vIn]; + + numTotalPoints = 0; + *intermCounts = 0; + + UInt32 faceCount = numFaces; + for( i = 0; i < faceCount; ++i ) + { + const UInt8 numInputPoints = *inputCounts; + if(ClipPolysByPlane(numInputPoints, inputPoints, planes[p], outputCounts, outputPoints, intermCounts, intermPoints)) + { + const UInt8 outputCount = *outputCounts++; + numTotalPoints+=outputCount; + outputPoints += outputCount; + } + else + { + if(0 == (--numFaces)) + break; + } + + inputCounts++; + inputPoints += numInputPoints; + } + vIn = 1-vIn; // anyone for ping-pong ? + + // add an extra face built from all the intersection points so it catches all the edges. + const UInt8 numIntermPoints = *intermCounts; + if(numIntermPoints && (p<5)) + { + numFaces++; + *outputCounts = numIntermPoints; + memcpy(outputPoints, intermPoints, numIntermPoints * sizeof(Vector3f)); + } + } + + if(numFaces) + { // output the clipped points + + Vector3f pt; + Vector3f ld = -lightDir; + points.reserve(numTotalPoints << 1); // worst case scenario + Vector3f* __restrict inputPoints=pData[vIn]; + UInt8* __restrict inputCounts=pDataCounts[vIn]; + + for( i = 0; i < numFaces; ++i ) + { + const UInt8 numPoints = *inputCounts++; + for(UInt8 k=0;k!=numPoints;k++) + { + const Vector3f& v = inputPoints[k]; + points.push_back(v); + if( IntersectLineAABB( pt, v, ld, sceneAABB ) ) + points.push_back( pt ); + } + inputPoints += numPoints; + } + } +} + +void CalculateBoundingSphereFromFrustumPoints(const Vector3f points[8], Vector3f& outCenter, float& outRadius) +{ + Vector3f spherePoints[4]; + spherePoints[0] = points[0]; + spherePoints[1] = points[3]; + spherePoints[2] = points[5]; + spherePoints[3] = points[7]; + + // Is bounding sphere at the far or near plane? + for( int plane = 1; plane >= 0; --plane ) + { + Vector3f pointA = spherePoints[plane*2]; + Vector3f pointB = spherePoints[plane*2 + 1]; + Vector3f center = (pointA + pointB) * 0.5f; + float radius2 = SqrMagnitude(pointA - center); + Vector3f pointC = spherePoints[(1-plane)*2]; + Vector3f pointD = spherePoints[(1-plane)*2 + 1]; + + // Check if all points are inside sphere + if( SqrMagnitude(pointC - center) <= radius2 && + SqrMagnitude(pointD - center) <= radius2 ) + { + outCenter = center; + outRadius = sqrt(radius2); + return; + } + } + + // Sphere touches all four frustum points + CalculateSphereFrom4Points(spherePoints, outCenter, outRadius); +} + +void CalculateSphereFrom4Points(const Vector3f points[4], Vector3f& outCenter, float& outRadius) +{ + Matrix4x4f mat; + + for( int i = 0; i < 4; ++i ) + { + mat.Get(i, 0) = points[i].x; + mat.Get(i, 1) = points[i].y; + mat.Get(i, 2) = points[i].z; + mat.Get(i, 3) = 1; + } + float m11 = mat.GetDeterminant(); + + for( int i = 0; i < 4; ++i ) + { + mat.Get(i, 0) = points[i].x*points[i].x + points[i].y*points[i].y + points[i].z*points[i].z; + mat.Get(i, 1) = points[i].y; + mat.Get(i, 2) = points[i].z; + mat.Get(i, 3) = 1; + } + float m12 = mat.GetDeterminant(); + + for( int i = 0; i < 4; ++i ) + { + mat.Get(i, 0) = points[i].x; + mat.Get(i, 1) = points[i].x*points[i].x + points[i].y*points[i].y + points[i].z*points[i].z; + mat.Get(i, 2) = points[i].z; + mat.Get(i, 3) = 1; + } + float m13 = mat.GetDeterminant(); + + for( int i = 0; i < 4; ++i ) + { + mat.Get(i, 0) = points[i].x; + mat.Get(i, 1) = points[i].y; + mat.Get(i, 2) = points[i].x*points[i].x + points[i].y*points[i].y + points[i].z*points[i].z; + mat.Get(i, 3) = 1; + } + float m14 = mat.GetDeterminant(); + + for( int i = 0; i < 4; ++i ) + { + mat.Get(i, 0) = points[i].x*points[i].x + points[i].y*points[i].y + points[i].z*points[i].z; + mat.Get(i, 1) = points[i].x; + mat.Get(i, 2) = points[i].y; + mat.Get(i, 3) = points[i].z; + } + float m15 = mat.GetDeterminant(); + + Vector3f c; + c.x = 0.5 * m12 / m11; + c.y = 0.5 * m13 / m11; + c.z = 0.5 * m14 / m11; + outRadius = sqrt(c.x*c.x + c.y*c.y + c.z*c.z - m15/m11); + outCenter = c; +} + +void CalculateSpotLightBounds (const float range, const float cotanHalfSpotAngle, const Matrix4x4f& lightMatrix, SpotLightBounds& outBounds) +{ + float sideLength = range / cotanHalfSpotAngle; + outBounds.points[0] = lightMatrix.GetPosition(); + outBounds.points[1] = lightMatrix.MultiplyPoint3( Vector3f(-sideLength,-sideLength, range) ); + outBounds.points[2] = lightMatrix.MultiplyPoint3( Vector3f( sideLength,-sideLength, range) ); + outBounds.points[3] = lightMatrix.MultiplyPoint3( Vector3f( sideLength, sideLength, range) ); + outBounds.points[4] = lightMatrix.MultiplyPoint3( Vector3f(-sideLength, sideLength, range) ); +} + +#if ENABLE_UNIT_TESTS + +#include "External/UnitTest++/src/UnitTest++.h" +#include "Runtime/Math/Random/rand.h" + +SUITE (BoundingUtilsTest) +{ +TEST(BoundingUtilsTest_CalculateSphereFrom4Points) +{ + Rand rnd(123); + for( int i = 0; i < 10; ++i ) + { + Vector3f points[4]; + for( int j = 0; j < 4; ++j ) + { + points[j].x = rnd.GetSignedFloat()*100.f; + points[j].y = rnd.GetSignedFloat()*100.f; + points[j].z = rnd.GetSignedFloat()*100.f; + } + Vector3f center; + float radius; + CalculateSphereFrom4Points(points, center, radius); + for( int j = 0; j < 4; ++j ) + { + float dist = Magnitude(points[j] - center); + float ratio = dist / radius; + // Radius may be large compared to input range, avoid comparing absolute distances + bool correct = (ratio > 0.999f) && (ratio < 1.001f); + CHECK(correct); + } + } +} +} + +#endif diff --git a/Runtime/Geometry/BoundingUtils.h b/Runtime/Geometry/BoundingUtils.h new file mode 100644 index 0000000..5338706 --- /dev/null +++ b/Runtime/Geometry/BoundingUtils.h @@ -0,0 +1,30 @@ +#pragma once + + +#include "Runtime/Math/Vector3.h" +#include <vector> + +class MinMaxAABB; +class Plane; +class Matrix4x4f; + +struct SpotLightBounds +{ + enum { kPointCount = 5 }; + Vector3f points[kPointCount]; +}; + +// 7-----6 far +// / /| +// 3-----2 | +// | 4 | 5 +// | |/ +// 0-----1 near + +void GetFrustumPoints( const Matrix4x4f& clipToWorld, Vector3f* frustum ); +void GetFrustumPortion( const Vector3f* frustum, float nearSplit, float farSplit, Vector3f* outPortion ); +void CalcHullBounds(const Vector3f* __restrict hullPoints, const UInt8* __restrict hullCounts, UInt8 hullFaces, const Plane& nearPlane, const Matrix4x4f& cameraWorldToClip, MinMaxAABB& aabb); +void CalculateFocusedLightHull( const Vector3f* frustum, const Vector3f& lightDir, const MinMaxAABB& sceneAABB, std::vector<Vector3f>& points ); +void CalculateBoundingSphereFromFrustumPoints(const Vector3f points[8], Vector3f& outCenter, float& outRadius); +void CalculateSphereFrom4Points(const Vector3f points[4], Vector3f& outCenter, float& outRadius); +void CalculateSpotLightBounds (const float range, const float cotanHalfSpotAngle, const Matrix4x4f& lightMatrix, SpotLightBounds& outBounds); diff --git a/Runtime/Geometry/BoundingVolumeConversion.h b/Runtime/Geometry/BoundingVolumeConversion.h new file mode 100644 index 0000000..059c163 --- /dev/null +++ b/Runtime/Geometry/BoundingVolumeConversion.h @@ -0,0 +1,32 @@ +#ifndef BOUNDINGVOLUMECONVERSION_H +#define BOUNDINGVOLUMECONVERSION_H + +#include "AABB.h" +#include "Sphere.h" + +inline void AABBToSphere (const AABB& a, Sphere& s) +{ + s.GetCenter () = a.GetCenter (); + s.GetRadius () = Magnitude (a.GetExtent ()); +} + +inline void MinMaxAABBToSphere (const MinMaxAABB& a, Sphere& s) +{ + AABB aabb = a; + AABBToSphere (aabb, s); +} + +inline void SphereToAABB (const Sphere& s, AABB& a) +{ + a.GetCenter () = s.GetCenter (); + a.GetExtent () = Vector3f (s.GetRadius (), s.GetRadius (), s.GetRadius ()); +} + +inline void SphereToMinMaxAABB (const Sphere& s, MinMaxAABB& minmax) +{ + AABB aabb; + SphereToAABB (s, aabb); + minmax = aabb; +} + +#endif diff --git a/Runtime/Geometry/ComputionalGeometry.cpp b/Runtime/Geometry/ComputionalGeometry.cpp new file mode 100644 index 0000000..5b979a5 --- /dev/null +++ b/Runtime/Geometry/ComputionalGeometry.cpp @@ -0,0 +1,304 @@ +#include "UnityPrefix.h" +#include "ComputionalGeometry.h" +#include "Runtime/Math/FloatConversion.h" +#include "Runtime/Math/Vector2.h" +#include "Plane.h" +#include <cmath> + +float SignedTriangleArea2D (const Vector2f& a, const Vector2f& b, const Vector2f& c) +{ + float i01 = (b.x - a.x) * (b.y + a.y); + float i12 = (c.x - b.x) * (c.y + b.y); + float i20 = (a.x - c.x) * (a.y + c.y); + + return (i01 + i12 + i20) * 0.5F; +} + +float SignedTriangleArea2D (const Vector3f* v) +{ + float i01 = (v[1].x - v[0].x) * (v[1].y + v[0].y); + float i12 = (v[2].x - v[1].x) * (v[2].y + v[1].y); + float i20 = (v[0].x - v[2].x) * (v[0].y + v[2].y); + + return (i01 + i12 + i20) * 0.5F; +} + +float TriangleArea3D (const Vector3f& a, const Vector3f& b, const Vector3f& c) +{ + return Magnitude (CalcRawNormalFromTriangle (a, b, c)) * 0.5F; +} + + +float CalculateProjectedBoxArea2D (const Vector3f* v) +{ + float i01 = (v[1].x - v[0].x) * (v[1].y + v[0].y); + float i13 = (v[3].x - v[1].x) * (v[3].y + v[1].y); + float i20 = (v[0].x - v[2].x) * (v[0].y + v[2].y); + float i37 = (v[7].x - v[3].x) * (v[7].y + v[3].y); + float i62 = (v[2].x - v[6].x) * (v[2].y + v[6].y); + float i23 = (v[3].x - v[2].x) * (v[3].y + v[2].y); + float i45 = (v[5].x - v[4].x) * (v[5].y + v[4].y); + float i57 = (v[7].x - v[5].x) * (v[7].y + v[5].y); + float i76 = (v[6].x - v[7].x) * (v[6].y + v[7].y); + float i64 = (v[4].x - v[6].x) * (v[4].y + v[6].y); + float i15 = (v[5].x - v[1].x) * (v[5].y + v[1].y); + float i40 = (v[0].x - v[4].x) * (v[0].y + v[4].y); + + float area + = Abs ( + i01 + + i13 + - i23 + + i20); + + area += Abs ( + i23 + + i37 + + i76 + + i62); + + area += Abs ( + i45 + + i57 + + i76 + + i64); + + area += Abs ( + i01 + + i15 + - i45 + + i40); + + area += Abs ( + i15 + + i57 + - i37 + - i13); + + area += Abs ( + - i40 + - i64 + + i62 + + i20); + + return area * 0.5F; +} + +float CalculateBoxAreaRadialFrustum ( + Vector3f* v, + float inFovy, + float inNear, + float inFar, + float inScreenHeight) +{ + UInt32 i; + double cotan, w, distance; + Vector3f projectedTriangle[8]; + + cotan = Deg2Rad (inFovy / 2.0F); + cotan = cos (cotan) / sin (cotan); + + for (i=0;i<8;i++) + { + projectedTriangle[i] = v[i]; + distance = Magnitude (projectedTriangle[i]); + + w = (inFar + inNear) / (inFar - inNear) * distance + + 2.0F * inNear * inFar / (inFar - inNear); + w = 1.0F / w * cotan * inScreenHeight / 2.0F; + + projectedTriangle[i].x *= w; + projectedTriangle[i].y *= w; + } + + return CalculateProjectedBoxArea2D (projectedTriangle); +} + +//float a = (inFar + inNear) / (inFar - inNear); +//float b = 2.0F * inNear * inFar / (inFar - inNear); +//float c = cotan * inScreenHeight / 2.0F; + +float CalculateBoxAreaRadialFrustum2 (Vector3f* v, float a, float b, float c) +{ + UInt32 i; + float w, distance; + Vector3f projectedTriangle[8]; + + for (i=0;i<8;i++) + { + projectedTriangle[i] = v[i]; + distance = Magnitude (projectedTriangle[i]); + + w = a * distance + b; + w = 1.0F / w * c; + + projectedTriangle[i].x *= w; + projectedTriangle[i].y *= w; + } + + return CalculateProjectedBoxArea2D (projectedTriangle); + +} + +float CalculateTriangleAreaRadialFrustum ( + Vector3f* v, + float inFovy, + float inNear, + float inFar, + float inScreenHeight) +{ + UInt32 i; + double cotan, w, distance; + Vector3f projectedTriangle[3]; + + cotan = Deg2Rad (inFovy / 2.0F); + cotan = cos (cotan) / sin (cotan); + + for (i=0;i<3;i++) + { + projectedTriangle[i] = v[i]; + distance = Magnitude (projectedTriangle[i]); + + w = (inFar + inNear) / (inFar - inNear) * distance + + 2.0F * inNear * inFar / (inFar - inNear); + w = 1.0F / w; + + projectedTriangle[i].x *= w * cotan * inScreenHeight / 2.0F; + projectedTriangle[i].y *= w * cotan * inScreenHeight / 2.0F; + } + + return Abs (SignedTriangleArea2D (projectedTriangle)); +} +//float a = (inFar + inNear) / (inFar - inNear); +//float b = 2.0F * inNear * inFar / (inFar - inNear); +//float c = cotan * inScreenHeight / 2.0F; + +float CalculateTriangleAreaRadialFrustum2 (Vector3f* v, float a, float b, float c) +{ + UInt32 i; + float w, distance; + Vector3f projectedTriangle[3]; + + for (i=0;i<3;i++) + { + projectedTriangle[i] = v[i]; + distance = Magnitude (projectedTriangle[i]); + + w = a * distance + b; + w = 1.0F / w * c; + + projectedTriangle[i].x *= w; + projectedTriangle[i].y *= w; + } + + return Abs (SignedTriangleArea2D (projectedTriangle)); +} + +float CalculateTriangleAreaRotationless (Vector3f* v, float inFovy, float inScreenWidth, float inScreenHeight, Vector3f& viewPoint) +{ + Vector3f normal = Normalize (Cross (v[0] - v[1], v[1] - v[2])); + float objectSpaceArea = Magnitude (Cross (v[1] - v[0], v[2] - v[0])) * 0.5F; + Vector3f centroid = (v[0] + v[1] + v[2]) / 3.0F; + Vector3f difference = viewPoint - centroid; + float distance = Magnitude (difference); + // We want 1.0 when it faces the tri directly -> 0 degrees + // We want 0.0 when the tri is invisible -> 90 degrees + // It doesn't matter if the tri faces away from the viewer or not. + float angle = Abs (Dot (normal, (difference / distance))); + + float screenSpaceArea = + objectSpaceArea + * sin (Deg2Rad (inFovy * 0.5F)) + * inScreenWidth + * inScreenHeight + * sin (Deg2Rad (inFovy * 0.5F)) + * angle + / distance; + + return screenSpaceArea; +} + +int ClipPolygonAgainstPlane(int vertexCount, const Vector3f *vertex, const Plane& plane, char *location, Vector3f *result) +{ + const float kBoundaryEpsilon = 1.0e-3F; + + enum + { + kPolygonInterior = 1, //## The point lies in the interior of the polygon. + kPolygonBoundary = 0, //## The point lies on or very near the boundary of the polygon. + kPolygonExterior = -1 //## The point lies outside the polygon. + }; + + int positive = 0; + int negative = 0; + + for (int a = 0; a < vertexCount; a++) + { + float d = plane.GetDistanceToPoint(vertex[a]); + if (d > kBoundaryEpsilon) + { + location[a] = kPolygonInterior; + positive++; + } + else + { + if (d < -kBoundaryEpsilon) + { + location[a] = kPolygonExterior; + negative++; + } + else + { + location[a] = kPolygonBoundary; + } + } + } + + if (negative == 0) + { + for (int a = 0; a < vertexCount; a++) + result[a] = vertex[a]; + return vertexCount; + } + else if (positive == 0) + { + return 0; + } + + int count = 0; + int previous = vertexCount - 1; + for (int index = 0; index < vertexCount; index++) + { + int loc = location[index]; + if (loc == kPolygonExterior) + { + if (location[previous] == kPolygonInterior) + { + const Vector3f& v1 = vertex[previous]; + const Vector3f& v2 = vertex[index]; + Vector3f dv = v2 - v1; + + float t = plane.GetDistanceToPoint(v2) / plane.GetDistanceToPoint(dv); + result[count++] = v2 - dv * t; + } + } + else + { + const Vector3f& v1 = vertex[index]; + if ((loc == kPolygonInterior) && (location[previous] == kPolygonExterior)) + { + const Vector3f& v2 = vertex[previous]; + Vector3f dv = v2 - v1; + + float t = plane.GetDistanceToPoint(v2) / plane.GetDistanceToPoint(dv); + result[count++] = v2 - dv * t; + } + + result[count++] = v1; + } + + previous = index; + } + + return count; +} diff --git a/Runtime/Geometry/ComputionalGeometry.h b/Runtime/Geometry/ComputionalGeometry.h new file mode 100644 index 0000000..42f5e47 --- /dev/null +++ b/Runtime/Geometry/ComputionalGeometry.h @@ -0,0 +1,19 @@ +#ifndef COMPUTIONALGEOMETRY_H +#define COMPUTIONALGEOMETRY_H + +class Vector2f; +class Vector3f; +class Plane; + +float SignedTriangleArea2D (const Vector2f& a, const Vector2f& b, const Vector2f& c); +float SignedTriangleArea2D (const Vector3f* v); +float TriangleArea3D (const Vector3f& a, const Vector3f& b, const Vector3f& c); +float CalculateProjectedBoxArea2D (const Vector3f* v); +float CalculateTriangleAreaRotationless (Vector3f* v, float inFovy, float inScreenWidth, float inScreenHeight, Vector3f& viewPoint); +float CalculateTriangleAreaRadialFrustum (Vector3f* v, float inFovy, float inNear, float inFar, float inScreenHeight); +float CalculateBoxAreaRadialFrustum (Vector3f* v, float inFovy, float inNear, float inFar, float inScreenHeight); +float CalculateBoxAreaRadialFrustum2 (Vector3f* v, float a, float b, float c); +float CalculateTriangleAreaRadialFrustum2 (Vector3f* v, float a, float b, float c); +int ClipPolygonAgainstPlane(int vertexCount, const Vector3f *vertex, const Plane& plane, char *location, Vector3f *result); + +#endif diff --git a/Runtime/Geometry/Intersection.cpp b/Runtime/Geometry/Intersection.cpp new file mode 100644 index 0000000..d1ae9b9 --- /dev/null +++ b/Runtime/Geometry/Intersection.cpp @@ -0,0 +1,940 @@ +#include "UnityPrefix.h" +#include "Intersection.h" +#include "Ray.h" +#include "Plane.h" +#include "Sphere.h" +#include "AABB.h" +#include "Runtime/Utilities/LogAssert.h" +#include "TriTriIntersect.h" +#include "Runtime/Utilities/BitUtility.h" +#include "Runtime/Camera/CullingParameters.h" + + +using namespace std; + +bool IntersectRayTriangle (const Ray& ray, const Vector3f& a, const Vector3f& b, const Vector3f& c) +{ + float t; + return IntersectRayTriangle (ray, a, b, c, &t); +} + +bool IntersectRayTriangle (const Ray& ray, const Vector3f& a, const Vector3f& b, const Vector3f& c, float* outT) +{ + const float kMinDet = 1e-6f; + + float t, u, v; + Vector3f edge1, edge2, tvec, pvec, qvec; + float det, inv_det; + + /* find vectors for two edges sharing vert0 */ + edge1 = b - a; + edge2 = c - a; + + /* begin calculating determinant - also used to calculate U parameter */ + pvec = Cross (ray.GetDirection (), edge2); + + /* if determinant is near zero, ray lies in plane of triangle */ + det = Dot (edge1, pvec); + + if (Abs (det) < kMinDet) + return false; + + inv_det = 1.0F / det; + + /* calculate distance from vert0 to ray origin */ + tvec = ray.GetOrigin () - a; + + /* calculate U parameter and test bounds */ + u = Dot (tvec, pvec) * inv_det; + if (u < 0.0F || u > 1.0F) + return false; + + /* prepare to test V parameter */ + qvec = Cross (tvec, edge1); + + /* calculate V parameter and test bounds */ + v = Dot (ray.GetDirection (), qvec) * inv_det; + if (v < 0.0F || u + v > 1.0F) + return false; + + t = Dot (edge2, qvec) * inv_det; + if (t < 0.0F) + return false; + *outT = t; + + return true; +} + +bool IntersectRaySphere (const Ray& ray, const Sphere& inSphere) +{ + Vector3f dif = inSphere.GetCenter () - ray.GetOrigin (); + float d = Dot (dif, ray.GetDirection ()); + float lSqr = Dot (dif, dif); + float rSqr = Sqr (inSphere.GetRadius ()); + + if (d < 0.0F && lSqr > rSqr) + return false; + + float mSqr = lSqr - Sqr (d); + + if (mSqr > rSqr) + return false; + else + return true; +} +/* +bool IntersectRaySphere (const Ray& ray, const Sphere& inSphere, float* t) +{ + AssertIf (t == NULL); + + Vector3f dif = inSphere.GetCenter () - ray.GetOrigin (); + float d = Dot (dif, ray.GetDirection ()); + float lSqr = Dot (dif, dif); + float rSqr = Sqr (inSphere.GetRadius ()); + + if (d < 0.0F && lSqr > rSqr) + return false; + + float mSqr = lSqr - Sqr (d); + + if (mSqr > rSqr) + return false; + + float q = sqrt (rSqr - mSqr); + + // ray.origin is inside the ray so a negative intersection will be returned + *t = d - q; + + return true; +} +*/ +bool IntersectRaySphere (const Ray& ray, const Sphere& inSphere, float* t0, float* t1) +{ + AssertIf (t0 == NULL); + AssertIf (t1 == NULL); + + Vector3f dif = inSphere.GetCenter () - ray.GetOrigin (); + float d = Dot (dif, ray.GetDirection ()); + float lSqr = Dot (dif, dif); + float rSqr = Sqr (inSphere.GetRadius ()); + + if (d < 0.0F && lSqr > rSqr) + return false; + + float mSqr = lSqr - Sqr (d); + + if (mSqr > rSqr) + return false; + + float q = sqrt (rSqr - mSqr); + + *t0 = d - q; + *t1 = d + q; + + return true; +} + +bool IntersectRayAABB (const Ray& ray, const AABB& inAABB) +{ + float t0, t1; + return IntersectRayAABB (ray, inAABB, &t0, &t1); +} +/* +bool IntersectRayAABB (const Ray& ray, const AABB& inAABB, float* outT) +{ + float tmin = -Vector3f::infinity; + float tmax = Vector3f::infinity; + + float t0, t1, f; + + Vector3f p = inAABB.GetCenter () - ray.GetOrigin (); + Vector3f extent = inAABB.GetExtent (); + long i; + for (i=0;i<3;i++) + { + // ray and plane are paralell so no valid intersection can be found + //if (Abs (ray.GetDirection ()[i]) > Vector3f::epsilon) + { + f = 1.0F / ray.GetDirection ()[i]; + t0 = (p[i] + extent[i]) * f; + t1 = (p[i] - extent[i]) * f; + // Ray leaves on Right, Top, Back Side + if (t0 < t1) + { + if (t0 > tmin) + tmin = t0; + + if (t1 < tmax) + tmax = t1; + + if (tmin > tmax) + return false; + + if (tmax < 0.0F) + return false; + } + // Ray leaves on Left, Bottom, Front Side + else + { + if (t1 > tmin) + tmin = t1; + + if (t0 < tmax) + tmax = t0; + + if (tmin > tmax) + return false; + + if (tmax < 0.0F) + return false; + } + } + } + + if (tmin > 0.0F) + *outT = tmin; + // ray starts inside the aabb + else + *outT = 0.0F; + + AssertIf (*outT < 0.0F); + return true; +} +*/ + +bool IntersectRayAABB (const Ray& ray, const AABB& inAABB, float* outT0) +{ + float t1; + return IntersectRayAABB (ray, inAABB, outT0, &t1); +} + +bool IntersectRayAABB (const Ray& ray, const AABB& inAABB, float* outT0, float* outT1) +{ + float tmin = -Vector3f::infinity; + float tmax = Vector3f::infinity; + + float t0, t1, f; + + Vector3f p = inAABB.GetCenter () - ray.GetOrigin (); + Vector3f extent = inAABB.GetExtent (); + long i; + for (i=0;i<3;i++) + { + // ray and plane are paralell so no valid intersection can be found + { + f = 1.0F / ray.GetDirection ()[i]; + t0 = (p[i] + extent[i]) * f; + t1 = (p[i] - extent[i]) * f; + // Ray leaves on Right, Top, Back Side + if (t0 < t1) + { + if (t0 > tmin) + tmin = t0; + + if (t1 < tmax) + tmax = t1; + + if (tmin > tmax) + return false; + + if (tmax < 0.0F) + return false; + } + // Ray leaves on Left, Bottom, Front Side + else + { + if (t1 > tmin) + tmin = t1; + + if (t0 < tmax) + tmax = t0; + + if (tmin > tmax) + return false; + + if (tmax < 0.0F) + return false; + } + } + } + + *outT0 = tmin; + *outT1 = tmax; + + return true; +} + + +bool IntersectSphereSphere (const Sphere& s0, const Sphere& s1) +{ + float sqrDist = SqrMagnitude (s0.GetCenter () - s1.GetCenter ()); + if (Sqr (s0.GetRadius () + s1.GetRadius ()) > sqrDist) + return true; + else + return false; +} + +bool IntersectSphereSphereInclusive (const Sphere& s0, const Sphere& s1) +{ + float sqrDist = SqrMagnitude (s0.GetCenter () - s1.GetCenter ()); + if (Sqr (s0.GetRadius () + s1.GetRadius ()) >= sqrDist) + return true; + else + return false; +} + +bool IntersectAABBAABB (const AABB& b0, const AABB& b1) +{ + const Vector3f dif = (b1.GetCenter () - b0.GetCenter ()); + + return Abs (dif.x) < b0.GetExtent (0) + b1.GetExtent (0) + && Abs (dif.y) < b0.GetExtent (1) + b1.GetExtent (1) + && Abs (dif.z) < b0.GetExtent (2) + b1.GetExtent (2); +} + +bool IntersectAABBAABBInclusive (const AABB& b0, const AABB& b1) +{ + const Vector3f dif = (b1.GetCenter () - b0.GetCenter ()); + + return Abs (dif.x) <= b0.GetExtent (0) + b1.GetExtent (0) + && Abs (dif.y) <= b0.GetExtent (1) + b1.GetExtent (1) + && Abs (dif.z) <= b0.GetExtent (2) + b1.GetExtent (2); +} + +bool IntersectAABBSphere (const AABB& aabb, const Sphere& s) +{ + return CalculateSqrDistance (s.GetCenter (), aabb) < Sqr (s.GetRadius ()); +} + +bool IntersectAABBSphereInclusive (const AABB& aabb, const Sphere& s) +{ + return CalculateSqrDistance (s.GetCenter (), aabb) <= Sqr (s.GetRadius ()); +} + +// possible optimization: Precalculate the Abs () of the plane. (3 fabs less per plane) +bool IntersectAABBFrustum ( + const AABB& a, + const Plane* p, + UInt32 inClipMask) +{ + const Vector3f& m = a.GetCenter ();// center of AABB + const Vector3f& extent = a.GetExtent ();// half-diagonal + UInt32 mk = 1; + + // loop while there are active planes.. + while (mk <= inClipMask) + { + // if clip plane is active... + if (inClipMask & mk) + { + const Vector3f& normal = p->GetNormal (); + float dist = p->GetDistanceToPoint (m); + float radius = Dot (extent, Abs (normal)); + + if (dist + radius < 0) return false; // behind clip plane + // if (dist - radius < 0) *outClipMask |= mk; // straddles clipplane + // else in front of clip plane-> leave outClipMask bit off + // float m = (p->a () * b->v[p->nx].x) + (p->b * b->v[p->ny].y) + (p->c * b->v[p->nz].z); + // if (m > -p->d ()) return OUTSIDE; + // float r = Dot (m, normal) + p->d (); + // float n = (extent.x * Abs(normal.x)) + (extent.y * Abs(normal.y)) + (extent.z * Abs(normal.z)); + // if (r + n < 0) return false; + } + mk += mk; + p++; // next plane + } + return true; // AABB intersects frustum +} +// Optimize like this: http://www.cg.tuwien.ac.at/studentwork/CESCG/CESCG-2002/DSykoraJJelinek/ + +bool IntersectAABBFrustumFull (const AABB& a, const Plane p[6]) +{ + return IntersectAABBPlaneBounds (a, p, 6); +} + +bool IntersectAABBPlaneBounds (const AABB& a, const Plane* p, const int planeCount) +{ + const Vector3f& m = a.GetCenter ();// center of AABB + const Vector3f& extent = a.GetExtent ();// half-diagonal + + for (int i = 0; i < planeCount; ++i, ++p) + { + const Vector3f& normal = p->GetNormal (); + float dist = p->GetDistanceToPoint (m); + float radius = Dot (extent, Abs (normal)); + if (dist + radius < 0) return false; // behind clip plane + } + return true; // AABB intersects space bounded by planes +} + +// Returns the shortest distance to planes if point is outside (positive float), +// and 0.0 if point is inside frustum planes +float PointDistanceToFrustum (const Vector4f& point, const Plane* p, const int planeCount) +{ + float maxDistanceNegative = -std::numeric_limits<float>::infinity(); + for (int i = 0; i < planeCount; ++i, ++p) + { + float dist = p->GetDistanceToPoint (point); + if ((dist<0.0f) && (dist > maxDistanceNegative)) + maxDistanceNegative = dist; + } + if (maxDistanceNegative != -std::numeric_limits<float>::infinity()) + return -maxDistanceNegative; + else + return 0.0f; +} + +bool IntersectTriTri (const Vector3f& a0, const Vector3f& b0, const Vector3f& c0, + const Vector3f& a1, const Vector3f& b1, const Vector3f& c1, + Vector3f* intersectionLine0, Vector3f* intersectionLine1, bool* coplanar) +{ + int coplanarInt; + bool ret; + ret = tri_tri_intersect_with_isectline (const_cast<Vector3f&> (a0).GetPtr (), const_cast<Vector3f&> (b0).GetPtr (), const_cast<Vector3f&> (c0).GetPtr (), + const_cast<Vector3f&> (a1).GetPtr (), const_cast<Vector3f&> (b1).GetPtr (), const_cast<Vector3f&> (c1).GetPtr (), + &coplanarInt, intersectionLine0->GetPtr (), intersectionLine1->GetPtr ()); + *coplanar = coplanarInt; + return ret; +} + +bool IntersectRayPlane (const Ray& ray, const Plane& plane, float* enter) +{ + AssertIf (enter == NULL); + float vdot = Dot (ray.GetDirection (), plane.GetNormal ()); + float ndot = -Dot (ray.GetOrigin (), plane.GetNormal ()) - plane.d (); + + // is line parallel to the plane? if so, even if the line is + // at the plane it is not considered as intersection because + // it would be impossible to determine the point of intersection + if ( CompareApproximately (vdot, 0.0F) ) + return false; + + // the resulting intersection is behind the origin of the ray + // if the result is negative ( enter < 0 ) + *enter = ndot / vdot; + + return *enter > 0.0F; +} + +bool IntersectSegmentPlane( const Vector3f& p1, const Vector3f& p2, const Plane& plane, Vector3f* result ) +{ + AssertIf( result == NULL ); + Vector3f vec = p2 - p1; + float vdot = Dot( vec, plane.GetNormal() ); + + // segment parallel to the plane + if( CompareApproximately(vdot, 0.0f) ) + return false; + + float ndot = -Dot( p1, plane.GetNormal() ) - plane.d(); + float u = ndot / vdot; + // intersection is out of segment + if( u < 0.0f || u > 1.0f ) + return false; + + *result = p1 + vec * u; + return true; +} + +/* +bool IntersectSphereTriangle (const Sphere& s, const Vector3f& vert0, const Vector3f& vert1, const Vector3f& vert2) +{ + const Vector3f& center = s.GetCenter (); + float radius = s.GetRadius (); + float radius2 = radius * radius; + // Early exit if one of the vertices is inside the sphere + Vector3f kDiff = vert2 - center; + float fC = SqrMagnitude (kDiff); + if(fC <= radius2) return true; + + kDiff = vert1 - center; + fC = SqrMagnitude (kDiff); + if(fC <= radius2) return true; + + kDiff = vert0 - center; + fC = SqrMagnitude (kDiff); + if(fC <= radius2) return true; + + // Else do the full distance test + Vector3f TriEdge0 = vert1 - vert0; + Vector3f TriEdge1 = vert2 - vert0; + + float fA00 = SqrMagnitude (TriEdge0); + float fA01 = Dot (TriEdge0, TriEdge1); + float fA11 = SqrMagnitude (TriEdge1); + float fB0 = Dot (kDiff, TriEdge0); + float fB1 = Dot (kDiff, TriEdge0); + + float fDet = Abs(fA00*fA11 - fA01*fA01); + float u = fA01*fB1-fA11*fB0; + float v = fA01*fB0-fA00*fB1; + float SqrDist; + + if(u + v <= fDet) + { + if(u < 0.0f) + { + if(v < 0.0f) // region 4 + { + if(fB0 < 0.0f) + { + if(-fB0>=fA00) { SqrDist = fA00+2.0f*fB0+fC; } + else { u = -fB0/fA00; SqrDist = fB0*u+fC; } + } + else + { + if(fB1>=0.0f) { SqrDist = fC; } + else if(-fB1>=fA11) { SqrDist = fA11+2.0f*fB1+fC; } + else { v = -fB1/fA11; SqrDist = fB1*v+fC; } + } + } + else // region 3 + { + if(fB1>=0.0f) { SqrDist = fC; } + else if(-fB1>=fA11) { SqrDist = fA11+2.0f*fB1+fC; } + else { v = -fB1/fA11; SqrDist = fB1*v+fC; } + } + } + else if(v < 0.0f) // region 5 + { + if(fB0>=0.0f) { SqrDist = fC; } + else if(-fB0>=fA00) { SqrDist = fA00+2.0f*fB0+fC; } + else { u = -fB0/fA00; SqrDist = fB0*u+fC; } + } + else // region 0 + { + // minimum at interior point + if(fDet==0.0f) + { + SqrDist = std::numeric_limits<float>::max (); + } + else + { + float fInvDet = 1.0f/fDet; + u *= fInvDet; + v *= fInvDet; + SqrDist = u*(fA00*u+fA01*v+2.0f*fB0) + v*(fA01*u+fA11*v+2.0f*fB1)+fC; + } + } + } + else + { + float fTmp0, fTmp1, fNumer, fDenom; + + if(u < 0.0f) // region 2 + { + fTmp0 = fA01 + fB0; + fTmp1 = fA11 + fB1; + if(fTmp1 > fTmp0) + { + fNumer = fTmp1 - fTmp0; + fDenom = fA00-2.0f*fA01+fA11; + if(fNumer >= fDenom) + { +// u = 1.0f; +// v = 0.0f; + SqrDist = fA00+2.0f*fB0+fC; + } + else + { + u = fNumer/fDenom; + v = 1.0f - u; + SqrDist = u*(fA00*u+fA01*v+2.0f*fB0) + v*(fA01*u+fA11*v+2.0f*fB1)+fC; + } + } + else + { +// u = 0.0f; + if(fTmp1 <= 0.0f) { SqrDist = fA11+2.0f*fB1+fC; } + else if(fB1 >= 0.0f) { SqrDist = fC; } + else { v = -fB1/fA11; SqrDist = fB1*v+fC; } + } + } + else if(v < 0.0f) // region 6 + { + fTmp0 = fA01 + fB1; + fTmp1 = fA00 + fB0; + if(fTmp1 > fTmp0) + { + fNumer = fTmp1 - fTmp0; + fDenom = fA00-2.0f*fA01+fA11; + if(fNumer >= fDenom) + { + SqrDist = fA11+2.0f*fB1+fC; + } + else + { + v = fNumer/fDenom; + u = 1.0f - v; + SqrDist = u*(fA00*u+fA01*v+2.0f*fB0) + v*(fA01*u+fA11*v+2.0f*fB1)+fC; + } + } + else + { + if(fTmp1 <= 0.0f) { SqrDist = fA00+2.0f*fB0+fC; } + else if(fB0 >= 0.0f) { SqrDist = fC; } + else { u = -fB0/fA00; SqrDist = fB0*u+fC; } + } + } + else // region 1 + { + fNumer = fA11 + fB1 - fA01 - fB0; + if(fNumer <= 0.0f) + { + SqrDist = fA11+2.0f*fB1+fC; + } + else + { + fDenom = fA00-2.0f*fA01+fA11; + if(fNumer >= fDenom) + { + SqrDist = fA00+2.0f*fB0+fC; + } + else + { + u = fNumer/fDenom; + v = 1.0f - u; + SqrDist = u*(fA00*u+fA01*v+2.0f*fB0) + v*(fA01*u+fA11*v+2.0f*fB1)+fC; + } + } + } + } + + return Abs (SqrDist) < radius2; +}*/ + +bool IntersectSphereTriangle (const Sphere& s, const Vector3f& vert0, const Vector3f& vert1, const Vector3f& vert2) +{ + const Vector3f& center = s.GetCenter (); + float radius = s.GetRadius (); + float radius2 = radius * radius; + Vector3f Diff; + + // Early exit if one of the vertices is inside the sphere + float sqrDiff; + Diff = vert1 - center; + sqrDiff = SqrMagnitude (Diff); + if(sqrDiff <= radius2) return true; + + Diff = vert2 - center; + sqrDiff = SqrMagnitude (Diff); + if(sqrDiff <= radius2) return true; + + Diff = vert0 - center; + sqrDiff = SqrMagnitude (Diff); + if(sqrDiff <= radius2) return true; + + // Else do the full distance test + Vector3f Edge0 = vert1 - vert0; + Vector3f Edge1 = vert2 - vert0; + + float A00 = Dot (Edge0, Edge0); + float A01 = Dot (Edge0, Edge1); + float A11 = Dot (Edge1, Edge1); + + float B0 = Dot (Diff, Edge0); + float B1 = Dot (Diff, Edge1); + + float C = Dot (Diff, Diff); + + float Det = Abs (A00 * A11 - A01 * A01); + float u = A01 * B1 - A11 * B0; + float v = A01 * B0 - A00 * B1; + + float DistSq; + if (u + v <= Det) + { + if(u < 0.0F) + { + if(v < 0.0F) + { + // region 4 + if(B0 < 0.0F) + { + if (-B0 >= A00) + DistSq = A00 + 2.0F * B0 + C; + else + { + u = -B0 / A00; + DistSq = B0 * u + C; + } + } + else{ + if(B1 >= 0.0F) + DistSq = C; + else if(-B1 >= A11) + DistSq = A11 + 2.0F * B1 + C; + else + { + v = -B1 / A11; + DistSq = B1 * v + C; + } + } + } + else + { // region 3 + if(B1 >= 0.0F) + DistSq = C; + else if(-B1 >= A11) + DistSq = A11 + 2.0F * B1 + C; + else + { + v = -B1 / A11; + DistSq = B1 * v + C; + } + } + } + else if(v < 0.0F) + { // region 5 + if (B0 >= 0.0F) + DistSq = C; + else if (-B0 >= A00) + DistSq = A00 + 2.0F * B0 + C; + else + { + u = -B0 / A00; + DistSq = B0 * u + C; + } + } + else + { // region 0 + // minimum at interior point + if (Det == 0.0F) + DistSq = std::numeric_limits<float>::max (); + else + { + float InvDet = 1.0F / Det; + u *= InvDet; + v *= InvDet; + DistSq = u * (A00 * u + A01 * v + 2.0F * B0) + v * (A01 * u + A11 * v + 2.0F * B1) + C; + } + } + } + else{ + double Tmp0, Tmp1, Numer, Denom; + + if(u < 0.0F) + { + // region 2 + Tmp0 = A01 + B0; + Tmp1 = A11 + B1; + if (Tmp1 > Tmp0){ + Numer = Tmp1 - Tmp0; + Denom = A00 - 2.0F * A01 + A11; + if (Numer >= Denom) + DistSq = A00 + 2.0F * B0 + C; + else + { + u = Numer / Denom; + v = 1.0 - u; + DistSq = u * (A00 * u + A01 * v + 2.0F * B0) + v * (A01 * u + A11 * v + 2.0F * B1) + C; + } + } + else + { + if(Tmp1 <= 0.0F) + DistSq = A11 + 2.0F * B1 + C; + else if(B1 >= 0.0) + DistSq = C; + else + { + v = -B1 / A11; + DistSq = B1 * v + C; + } + } + } + else if(v < 0.0) + { // region 6 + Tmp0 = A01 + B1; + Tmp1 = A00 + B0; + if (Tmp1 > Tmp0) + { + Numer = Tmp1 - Tmp0; + Denom = A00 - 2.0F * A01 + A11; + if (Numer >= Denom) + DistSq = A11 + 2.0 * B1 + C; + else + { + v = Numer / Denom; + u = 1.0F - v; + DistSq = u * (A00 * u + A01 * v + 2.0F * B0) + v * (A01 * u + A11 * v + 2.0F * B1) + C; + } + } + else + { + if (Tmp1 <= 0.0F) + DistSq = A00 + 2.0F * B0 + C; + else if(B0 >= 0.0F) + DistSq = C; + else + { + u = -B0 / A00; + DistSq = B0 * u + C; + } + } + } + else + { + // region 1 + Numer = A11 + B1 - A01 - B0; + if (Numer <= 0.0F) + DistSq = A11 + 2.0F * B1 + C; + else + { + Denom = A00 - 2.0F * A01 + A11; + if (Numer >= Denom) + DistSq = A00 + 2.0F * B0 + C; + else + { + u = Numer / Denom; + v = 1.0F - u; + DistSq = u * (A00 * u + A01 * v + 2.0F * B0) + v * (A01 * u + A11 * v + 2.0F * B1) + C; + } + } + } + } + + return Abs (DistSq) <= radius2; +} + +bool TestPlanesAABB(const Plane* planes, const int planeCount, const AABB& bounds) +{ + UInt32 planeMask = 0; + if (planeCount == 6) + planeMask = 63; + else + { + for (int i = 0; i < planeCount; ++i) + planeMask |= 1 << i; + } + + return IntersectAABBFrustum (bounds, planes, planeMask); +} + +// -------------------------------------------------------------------------- + +#if ENABLE_UNIT_TESTS + +#include "External/UnitTest++/src/UnitTest++.h" +#include "Runtime/Math/Random/Random.h" +#include "Runtime/Profiler/TimeHelper.h" + +void GenerateUnitFrustumPlanes (Plane* p) +{ + p[0].SetABCD(-1.0, 0.0, 0.0,-1.0); + p[1].SetABCD( 1.0, 0.0, 0.0, 1.0); + p[2].SetABCD( 0.0,-1.0, 0.0,-1.0); + p[3].SetABCD( 0.0, 1.0, 0.0, 1.0); + p[4].SetABCD( 0.0, 0.0,-1.0,-1.0); + p[5].SetABCD( 0.0, 0.0, 1.0, 1.0); +} + +float PointDistanceToFrustumRef (const Vector3f& point, const Plane* p, const int planeCount) +{ + DebugAssert(planeCount <= 6); + + float maxDistanceNegative = -std::numeric_limits<float>::infinity(); + float distances[6]; + + // Point distances to frustum planes + for (int i=0; i<planeCount; i++) + distances[i] = p[i].GetDistanceToPoint (point); + + // Replace positive distances with negative infinity. This simplifies the shortest negative distance search to maximum operators + for (int i=0; i<planeCount; i++) + distances[i] = (distances[i] > 0.0f) ? -std::numeric_limits<float>::infinity() : distances[i]; + + // Find the shortest negative distance from the distance values + for (int i=0; i<planeCount; i++) + maxDistanceNegative = (distances[i] > maxDistanceNegative) ? distances[i] : maxDistanceNegative; + + // If maxNegativeDistance is negative infinity, all distance were positive and the point is inside the frustum. In that case, return 0.0, otherwise return the shortest distance (abs value) + return (maxDistanceNegative == -std::numeric_limits<float>::infinity()) ? 0.0f : -maxDistanceNegative; +} + +void TestComparePerformancePointDistanceToFrustum (); + +SUITE (IntersectionTests) +{ + TEST(PointDistanceToFrustum) + { + Plane unitFrustumPlanes[6]; + GenerateUnitFrustumPlanes(unitFrustumPlanes); + + Rand r (1); // fixed seed produces fixed random series + + for (int i=0; i<1000; i++) + { + // Random coordinates (2x unit frustum volume) + float x = (r.GetFloat() - 0.5f) * 3.0f; + float y = (r.GetFloat() - 0.5f) * 3.0f; + float z = (r.GetFloat() - 0.5f) * 3.0f; + Vector3f point3f(x,y,z); + Vector4f point4f(x,y,z,0.0f); + + float distanceRef = PointDistanceToFrustumRef(point3f, unitFrustumPlanes, 6); + float distance = PointDistanceToFrustum(point4f, unitFrustumPlanes, 6); + + if (distanceRef > 0.0F) + { + CHECK_EQUAL(distance, distanceRef); + } + else + { + CHECK(distance <= 0.0F); + } + } + } +} // SUITE + +static void TestPerformancePointDistanceToFrustum () +{ + Plane unitFrustumPlanes[6]; + GenerateUnitFrustumPlanes(unitFrustumPlanes); + + Rand r (1); // fixed seed produces fixed random series + float x = (r.GetFloat() - 0.5f) * 2.0f; + float y = (r.GetFloat() - 0.5f) * 2.0f; + float z = (r.GetFloat() - 0.5f) * 2.0f; + Vector4f pointFloat4(x,y,z,0.0f); + float distanceSum = 0.0f; + + ABSOLUTE_TIME time = START_TIME; + + for (int i=0; i<1000000; i++) + { + distanceSum += PointDistanceToFrustum(pointFloat4, unitFrustumPlanes, 6); + } + + printf_console("\n\nTime Impl: %f ms %f\n\n", GetElapsedTimeInSeconds (time) * 1000.0F, distanceSum); +} + +static void TestPerformancePointDistanceToFrustumRef () +{ + Plane unitFrustumPlanes[6]; + GenerateUnitFrustumPlanes(unitFrustumPlanes); + + Rand r (1); // fixed seed produces fixed random series + float x = (r.GetFloat() - 0.5f) * 2.0f; + float y = (r.GetFloat() - 0.5f) * 2.0f; + float z = (r.GetFloat() - 0.5f) * 2.0f; + Vector3f point3f(x,y,z); + float distanceSum = 0.0f; + + ABSOLUTE_TIME time = START_TIME; + + for (int i=0; i<1000000; i++) + { + distanceSum += PointDistanceToFrustumRef(point3f, unitFrustumPlanes, 6); + } + + printf_console("\n\nTime REF: %f (%f)ms\n\n", GetElapsedTimeInSeconds (time) * 1000.0F, distanceSum); +} + +void TestComparePerformancePointDistanceToFrustum () +{ + TestPerformancePointDistanceToFrustum (); + TestPerformancePointDistanceToFrustumRef (); +} + +#endif // ENABLE_UNIT_TESTS
\ No newline at end of file diff --git a/Runtime/Geometry/Intersection.h b/Runtime/Geometry/Intersection.h new file mode 100644 index 0000000..4c96099 --- /dev/null +++ b/Runtime/Geometry/Intersection.h @@ -0,0 +1,95 @@ +#ifndef INTERSECTION_H +#define INTERSECTION_H + +#include "Runtime/Math/Simd/math.h" + +class Ray; +class OptimizedRay; +class Sphere; +class AABB; +class Plane; + +class Vector3f; +class Vector4f; + +// Intersects a Ray with a triangle. +bool IntersectRayTriangle (const Ray& ray, const Vector3f& a, const Vector3f& b, const Vector3f& c); +// t is to be non-Null and returns the first intersection point of the ray (ray.o + t * ray.dir) +bool IntersectRayTriangle (const Ray& ray, const Vector3f& a, const Vector3f& b, const Vector3f& c, float* t); + +// Intersects a ray with a volume. +// Returns true if the ray stats inside the volume or in front of the volume +bool IntersectRaySphere (const Ray& ray, const Sphere& inSphere); +bool IntersectRayAABB (const Ray& ray, const AABB& inAABB); + +// Intersects a ray with a volume. +// Returns true if the ray stats inside the volume or in front of the volume +// t0 is the first, t1 the second intersection. Both have to be non-NULL. +// (t1 is always positive, t0 is negative if the ray starts inside the volume) +bool IntersectRayAABB (const Ray& ray, const AABB& inAABB, float* t0, float* t1); +bool IntersectRayAABB (const Ray& ray, const AABB& inAABB, float* t0); +bool IntersectRaySphere (const Ray& ray, const Sphere& inSphere, float* t0, float* t1); + +// Do these volumes intersect each other? +bool IntersectSphereSphere (const Sphere& s0, const Sphere& s1); +bool IntersectAABBAABB (const AABB& s0, const AABB& s1); +bool IntersectAABBSphere (const AABB& s0, const Sphere& s1); + +// Do these volumes intersect or touch each other? +bool IntersectSphereSphereInclusive (const Sphere& s0, const Sphere& s1); +bool EXPORT_COREMODULE IntersectAABBAABBInclusive (const AABB& s0, const AABB& s1); +bool IntersectAABBSphereInclusive (const AABB& s0, const Sphere& s1); + +// Tests if the aabb is inside any of the planes enabled by inClipMask +// The bitmask tells which planes have to be tested. (For 6 planes the bitmask is 63) +bool IntersectAABBFrustum (const AABB& a, const Plane* p, UInt32 inClipMask); +bool IntersectAABBFrustumFull (const AABB& a, const Plane p[6]); +bool IntersectAABBPlaneBounds (const AABB& a, const Plane* p, const int planeCount); + +float PointDistanceToFrustum (const Vector4f& point, const Plane* p, const int planeCount); + +bool IntersectTriTri (const Vector3f& a0, const Vector3f& b0, const Vector3f& c0, + const Vector3f& a1, const Vector3f& b1, const Vector3f& c1, + Vector3f* intersectionLine0, Vector3f* intersectionLine1, bool* coplanar); + +// Intersects a ray with a plane (The ray can hit the plane from front and behind) +// On return enter is the rays parameter where the intersection occurred. +bool IntersectRayPlane (const Ray& ray, const Plane& plane, float* enter); + + +// Intersects a line segment with a plane (can hit the plane from front and behind) +// Fill result point if intersected. +bool IntersectSegmentPlane( const Vector3f& p1, const Vector3f& p2, const Plane& plane, Vector3f* result ); + + +// Returns true if the triangle touches or is inside the triangle (a, b, c) +bool IntersectSphereTriangle (const Sphere& s, const Vector3f& a, const Vector3f& b, const Vector3f& c); + +/// Returns true if the bounding box is inside the planes or intersects any of the planes. +bool TestPlanesAABB(const Plane* planes, const int planeCount, const AABB& bounds); + +/// Projects point on a line. +template <typename T> +T ProjectPointLine(const T& point, const T& lineStart, const T& lineEnd) +{ + T relativePoint = point - lineStart; + T lineDirection = lineEnd - lineStart; + float length = Magnitude(lineDirection); + T normalizedLineDirection = lineDirection; + if (length > T::epsilon) + normalizedLineDirection /= length; + + float dot = Dot(normalizedLineDirection, relativePoint); + dot = clamp(dot, 0.0F, length); + + return lineStart + normalizedLineDirection * dot; +} + +/// Returns the distance to a line from a point. +template <typename T> +float DistancePointLine(const T& point, const T& lineStart, const T& lineEnd) +{ + return Magnitude(ProjectPointLine<T>(point, lineStart, lineEnd) - point); +} + +#endif diff --git a/Runtime/Geometry/IntersectionTests.cpp b/Runtime/Geometry/IntersectionTests.cpp new file mode 100644 index 0000000..1e384fc --- /dev/null +++ b/Runtime/Geometry/IntersectionTests.cpp @@ -0,0 +1,288 @@ +#include "UnityPrefix.h" +#include "Configuration/UnityConfigure.h" + +#if ENABLE_UNIT_TESTS + +#include "External/UnitTest++/src/UnitTest++.h" +#include "Runtime/Geometry/Intersection.h" +#include "Runtime/Geometry/AABB.h" +#include "Runtime/Geometry/Sphere.h" +#include "Runtime/Geometry/Ray.h" + + +SUITE (IntersectionTests) +{ + + +// AABB Ray (ray inside aabb) +TEST (RayInsideAABB) +{ + AABB aabb; + Ray ray; + float t0, t1; + bool result; + + aabb.GetCenter () = Vector3f (5.0F, 10.0F, 20.0F); + aabb.GetExtent () = Vector3f (5.0F, 10.0F, 20.0F); + + ray.SetOrigin (aabb.GetCenter ()); + ray.SetDirection (Vector3f::zAxis); + + CHECK (IntersectRayAABB (ray, aabb)); + + result = IntersectRayAABB (ray, aabb, &t0, &t1); + CHECK (result); + CHECK_CLOSE (t0, -20.0F, 0.000001F); + CHECK_CLOSE (t1, 20.0F, 0.000001F); + + ray.SetDirection (-Vector3f::zAxis); + result = IntersectRayAABB (ray, aabb, &t0, &t1); + CHECK (result); + CHECK_CLOSE (t0, -20.0F, 0.000001F); + CHECK_CLOSE (t1, 20.0F, 0.000001F); +} + + +// AABB Ray (ray doesn't hit) +TEST (RayOutsideAABB) +{ + AABB aabb; + Ray ray; + float t0, t1; + + aabb.GetCenter () = Vector3f (5.0F, 10.0F, 20.0F); + aabb.GetExtent () = Vector3f (5.0F, 10.0F, 20.0F); + + ray.SetOrigin (aabb.GetCenter () + Vector3f (5.0F, 10.0F, 20.01F)); + ray.SetDirection (Vector3f::zAxis); + + CHECK (!IntersectRayAABB (ray, aabb)); + CHECK (!IntersectRayAABB (ray, aabb, &t0, &t1)); +} + + +// AABB Ray (ray hits frontal) +TEST (RayHitsAABBFrontal) +{ + AABB aabb; + Ray ray; + float t0, t1; + bool result; + + aabb.GetCenter () = Vector3f (5.0F, 10.0F, 20.0F); + aabb.GetExtent () = Vector3f (5.0F, 10.0F, 20.0F); + + ray.SetOrigin (Vector3f (5.0F, 10.0F, 60.0F)); + ray.SetDirection (-Vector3f::zAxis); + + CHECK (IntersectRayAABB (ray, aabb)); + + result = IntersectRayAABB (ray, aabb, &t0, &t1); + CHECK (result); + CHECK_CLOSE (t0, 20.0F, 0.000001F); + CHECK_CLOSE (t1, 60.0F, 0.000001F); +} + + +// AABB Ray (ray hits backward) +TEST (RayHitsAABBBackward) +{ + AABB aabb; + Ray ray; + float t0, t1; + + aabb.GetCenter () = Vector3f (5.0F, 10.0F, 20.0F); + aabb.GetExtent () = Vector3f (5.0F, 10.0F, 20.0F); + + ray.SetOrigin (Vector3f (5.0F, 10.0F, 60.0F)); + ray.SetDirection (Vector3f (0.0F, 0.0F, 1.0F)); + + CHECK (!IntersectRayAABB (ray, aabb)); + CHECK (!IntersectRayAABB (ray, aabb, &t0, &t1)); +} + + +// Sphere Ray (ray inside) (sphere.origin in front of ray.origin) +TEST (SphereToRay1) +{ + Sphere sphere; + Ray ray; + float t0, t1; + bool result; + + sphere.GetCenter () = Vector3f (5.0F, 10.0F, 20.0F); + sphere.GetRadius () = 10.0F; + + ray.SetOrigin (Vector3f (5.0F, 10.0F, 25.0F)); + ray.SetDirection (Vector3f (0.0F, 0.0F, 1.0F)); + + CHECK (IntersectRaySphere (ray, sphere)); + + result = IntersectRaySphere (ray, sphere, &t0, &t1); + CHECK (result); + CHECK_CLOSE (t0, -15.0F, 0.000001F); + CHECK_CLOSE (t1, 5.0F, 0.000001F); +} + + +// Sphere Ray (ray inside) (sphere.origin in front of ray.origin) +TEST (SphereRay2) +{ + Sphere sphere; + Ray ray; + float t0, t1; + bool result; + + sphere.GetCenter () = Vector3f (5.0F, 10.0F, 20.0F); + sphere.GetRadius () = 10.0F; + + ray.SetOrigin (Vector3f (5.0F, 10.0F, 25.0F)); + ray.SetDirection (Vector3f (0.0F, 0.0F, -1.0F)); + + CHECK (IntersectRaySphere (ray, sphere)); + + result = IntersectRaySphere (ray, sphere, &t0, &t1); + CHECK (result); + CHECK_CLOSE (t0, -5.0F, 0.000001F); + CHECK_CLOSE (t1, 15.0F, 0.000001F); +} + + +// Sphere Ray (hits sphere frontal) +TEST (RayHitsSphereFrontal) +{ + Sphere sphere; + Ray ray; + float t0, t1; + bool result; + + sphere.GetCenter () = Vector3f (5.0F, 10.0F, 20.0F); + sphere.GetRadius () = 10.0F; + + ray.SetOrigin (Vector3f (5.0F, 10.0F, 0.0F)); + ray.SetDirection (Vector3f (0.0F, 0.0F, 1.0F)); + + CHECK (IntersectRaySphere (ray, sphere)); + + result = IntersectRaySphere (ray, sphere, &t0, &t1); + CHECK (result); + CHECK_CLOSE (t0, 10.0F, 0.000001F); + CHECK_CLOSE (t1, 30.0F, 0.000001F); +} + + +// Sphere Ray (hits sphere backwards) +TEST (RayHitsSphereBackwards) +{ + Sphere sphere; + Ray ray; + float t0, t1; + bool result; + + sphere.GetCenter () = Vector3f (5.0F, 10.0F, 20.0F); + sphere.GetRadius () = 10.0F; + + ray.SetOrigin (Vector3f (5.0F, 10.0F, 40.0F)); + ray.SetDirection (Vector3f (0.0F, 0.0F, 1.0F)); + + CHECK (!IntersectRaySphere (ray, sphere)); + + result = IntersectRaySphere (ray, sphere, &t0, &t1); + CHECK (!result); + ErrorIf (result != false); +} + + +// Sphere Ray (misses sphere) +TEST (RayMissesSphere) +{ + Sphere sphere; + Ray ray; + float t0, t1; + + sphere.GetCenter () = Vector3f (5.0F, 10.0F, 20.0F); + sphere.GetRadius () = 10.0F; + + ray.SetOrigin (Vector3f (5.0F, 10.0F, 30.01F)); + ray.SetDirection (Vector3f (0.0F, 1.0F, 0.0F)); + + CHECK (!IntersectRaySphere (ray, sphere)); + CHECK (!IntersectRaySphere (ray, sphere, &t0, &t1)); +} + + +// Triangle Triangle (Not intersecting) +TEST (TriangleTriangleNotIntersecting) +{ + Vector3f + a1 (0, 0, 0), + a2 (1, 1, 0), + a3 (2, 0, 0), + b1 (0, 0, 1), + b2 (1, 1, 1), + b3 (2, 0, 1); + Vector3f r1, r2; + bool coplanar; + + CHECK (!IntersectTriTri (a1, a2, a3, b1, b2, b3, &r1, &r2, &coplanar)); +} + + +// Triangle Triangle (intersecting) +TEST (TriangleTriangleIntersecting) +{ + Vector3f + a1 (0, 2, 5), + a2 (2, 2, 0), + a3 (0, 2, 0), + b1 (0, 0, 0), + b2 (0, 5, 0), + b3 (0, 5, 3); + Vector3f r1, r2; + bool coplanar; + + CHECK (IntersectTriTri (a1, a2, a3, b1, b2, b3, &r1, &r2, &coplanar)); + CHECK (CompareApproximately (r1, Vector3f (0, 2, 0))); + CHECK (CompareApproximately (r2, Vector3f (0, 2, 1.2f))); + CHECK_EQUAL (false, coplanar); +} + + +// Triangle Triangle (coplanaer) +TEST (TriangleTriangleCoplanar) +{ + Vector3f + a1 (0, 8, 0), + a2 (0, 4, 0), + a3 (5, 4, 0), + b1 (0, 5, 0), + b2 (5, 0, 0), + b3 (0, 0, 0); + Vector3f r1, r2; + bool coplanar; + + CHECK (IntersectTriTri (a1, a2, a3, b1, b2, b3, &r1, &r2, &coplanar)); + CHECK (coplanar); +} + + +TEST (Misc) +{ + CHECK (IntersectSphereTriangle ( + Sphere (Vector3f (0.3F, 0.3F, 0.1F), .2F), + Vector3f (0.0F, 0.0F, 0.0F), + Vector3f (0.0F, 1.0F, 0.0F), + Vector3f (1.0F, 0.0F, 0.0F))); + CHECK (IntersectSphereTriangle ( + Sphere (Vector3f (0.3F, 0.3F, 0.0F), .2F), + Vector3f (0.0F, 0.0F, 0.0F), + Vector3f (0.0F, 1.0F, 0.0F), + Vector3f (1.0F, 0.0F, 0.0F))); + CHECK (!CompareApproximately (0.01f, 0.0)); +} + + +} + + +#endif diff --git a/Runtime/Geometry/Plane.h b/Runtime/Geometry/Plane.h new file mode 100644 index 0000000..ca58ff2 --- /dev/null +++ b/Runtime/Geometry/Plane.h @@ -0,0 +1,169 @@ +#ifndef PLANE_H +#define PLANE_H + +#include "Runtime/Math/Vector3.h" +#include "Runtime/Math/Vector4.h" +#include "Runtime/Serialize/SerializeUtility.h" + +// Calculates the unnormalized normal from 3 points +Vector3f CalcRawNormalFromTriangle (const Vector3f& a, const Vector3f& b, const Vector3f& c); + +enum{ + kPlaneFrustumLeft, + kPlaneFrustumRight, + kPlaneFrustumBottom, + kPlaneFrustumTop, + kPlaneFrustumNear, + kPlaneFrustumFar, + kPlaneFrustumNum, +}; + +class Plane +{ +public: + + Vector3f normal; + float distance; + + DECLARE_SERIALIZE_OPTIMIZE_TRANSFER (Plane) + + const float& a ()const { return normal.x; } + const float& b ()const { return normal.y; } + const float& c ()const { return normal.z; } + + const float& d () const { return distance; } + float& d () { return distance; } + + const Vector3f& GetNormal ()const{ return normal; } + + Plane () { } + Plane (float a, float b, float c, float d) { normal.x = a; normal.y = b; normal.z = c; distance = d; } + + Plane& operator *= (float scale); + bool operator == (const Plane& p)const { return normal == p.normal && distance == p.distance; } + bool operator != (const Plane& p)const { return normal != p.normal || distance != p.distance; } + + void SetInvalid () { normal = Vector3f::zero; distance = 0.0F; } + + // Just sets the coefficients. Does NOT normalize! + void SetABCD (const float a, const float b, const float c, const float d); + + void Set3Points (const Vector3f& a, const Vector3f& b, const Vector3f& c); + void Set3PointsUnnormalized (const Vector3f& a, const Vector3f& b, const Vector3f& c); + + void SetNormalAndPosition (const Vector3f& inNormal, const Vector3f& inPoint); + + float GetDistanceToPoint (const Vector3f& inPt) const; + float GetDistanceToPoint (const Vector4f& inPt) const; + bool GetSide (const Vector3f& inPt) const; + bool SameSide (const Vector3f& inPt0, const Vector3f& inPt1); + + void NormalizeRobust (); + void NormalizeUnsafe (); +}; + +template<class TransferFunction> +inline void Plane::Transfer (TransferFunction& transfer) +{ + TRANSFER (normal); + TRANSFER (distance); +} + +inline float Plane::GetDistanceToPoint (const Vector3f& inPt)const +{ + DebugAssert (IsNormalized (normal)); + return Dot (normal, inPt) + distance; +} + +// inPt w component is ignored from distance computations +inline float Plane::GetDistanceToPoint (const Vector4f& inPt)const +{ + DebugAssert (IsNormalized (normal)); + //Dot3 (normal, inPt) + distance; + return normal.x * inPt.x + normal.y * inPt.y + normal.z * inPt.z + distance; +} + +// Returns true if we are on the front side (same as: GetDistanceToPoint () > 0.0) +inline bool Plane::GetSide (const Vector3f& inPt) const +{ + return Dot (normal, inPt) + distance > 0.0F; +} + + +// Calculates the normal from 3 points unnormalized +inline Vector3f CalcRawNormalFromTriangle (const Vector3f& a, const Vector3f& b, const Vector3f& c) +{ + return Cross (b - a, c - a); +} + +inline Vector3f RobustNormalFromTriangle (const Vector3f& v0, const Vector3f& v1, const Vector3f& v2) +{ + Vector3f normal = CalcRawNormalFromTriangle (v0, v1, v2); + return NormalizeRobust (normal); +} + +inline void Plane::SetABCD (float a, float b, float c, float d) +{ + normal.Set(a, b, c); + distance = d; +} + +inline void Plane::Set3Points (const Vector3f& a, const Vector3f& b, const Vector3f& c) +{ + normal = CalcRawNormalFromTriangle (a, b, c); + normal = ::Normalize (normal); + distance = -Dot (normal, a); + AssertIf (!IsNormalized (normal)); +} + +inline void Plane::Set3PointsUnnormalized (const Vector3f& a, const Vector3f& b, const Vector3f& c) +{ + normal = CalcRawNormalFromTriangle (a, b, c); + distance = -Dot (normal, a); +} + +inline void Plane::SetNormalAndPosition (const Vector3f& inNormal, const Vector3f& inPoint) +{ + normal = inNormal; + AssertIf (!IsNormalized (normal, 0.001f)); + distance = -Dot (inNormal, inPoint); +} + +inline bool Plane::SameSide (const Vector3f& inPt0, const Vector3f& inPt1) +{ + float d0 = GetDistanceToPoint(inPt0); + float d1 = GetDistanceToPoint(inPt1); + if (d0 > 0.0f && d1 > 0.0f) + return true; + else if (d0 <= 0.0f && d1 <= 0.0f) + return true; + else + return false; +} + +inline Plane& Plane::operator *= (float scale) +{ + normal *= scale; + distance *= scale; + return *this; +} + +inline void Plane::NormalizeUnsafe () +{ + float invMag = 1.0f/Magnitude(normal); + normal *= invMag; + distance *= invMag; +} + +// It uses NormalizeRobust(), so it handles zero and extremely small vectors, +// but can be slow. Another option would be to use plain normalize, but +// always remember to check for division by zero with zero vectors. +inline void Plane::NormalizeRobust () +{ + float invMag; + normal = ::NormalizeRobust(normal, invMag); + distance *= invMag; +} + + +#endif diff --git a/Runtime/Geometry/Ray.cpp b/Runtime/Geometry/Ray.cpp new file mode 100644 index 0000000..0dac305 --- /dev/null +++ b/Runtime/Geometry/Ray.cpp @@ -0,0 +1,17 @@ +#include "UnityPrefix.h" +#include "Ray.h" + +float Ray::SqrDistToPoint(const Vector3f& P) const { +//ÊÊÊ Vector v = L.P1 - L.P0; + Vector3f v = m_Direction; +// Vector3f w = P - L.P0; + Vector3f w = P - m_Origin; + + float c1 = Dot(w,v); + float c2 = Dot(v,v); + float b = c1 / c2; + +//ÊÊÊ Point Pb = L.P0 + b * v; + Vector3f Pb = GetPoint (b); + return SqrMagnitude(P - Pb); +} diff --git a/Runtime/Geometry/Ray.h b/Runtime/Geometry/Ray.h new file mode 100644 index 0000000..515126f --- /dev/null +++ b/Runtime/Geometry/Ray.h @@ -0,0 +1,32 @@ +#ifndef RAY_H +#define RAY_H + +#include "Runtime/Math/FloatConversion.h" +#include "Runtime/Math/Vector3.h" + +class Ray +{ +#if UNITY_FLASH //flash needs to be able to set these fields +public: +#endif + Vector3f m_Origin; + Vector3f m_Direction; // Direction is always normalized + +public: + Ray () {} + Ray (const Vector3f& orig, const Vector3f& dir) { m_Origin = orig; SetDirection (dir); } + + const Vector3f& GetDirection ()const { return m_Direction; } + // Direction has to be normalized + void SetDirection (const Vector3f& dir) { AssertIf (!IsNormalized (dir)); m_Direction = dir; } + void SetApproxDirection (const Vector3f& dir) { m_Direction = NormalizeFast (dir); } + void SetOrigin (const Vector3f& origin) { m_Origin = origin; } + + const Vector3f& GetOrigin ()const { return m_Origin; } + Vector3f GetPoint (float t) const { return m_Origin + t * m_Direction; } + + float SqrDistToPoint (const Vector3f &v) const; +}; + + +#endif diff --git a/Runtime/Geometry/Ray2D.h b/Runtime/Geometry/Ray2D.h new file mode 100644 index 0000000..7adff8e --- /dev/null +++ b/Runtime/Geometry/Ray2D.h @@ -0,0 +1,32 @@ +#ifndef RAY2D_H +#define RAY2D_H + +#include "Runtime/Math/Vector2.h" + + +// -------------------------------------------------------------------------- + + +class Ray2D +{ +#if UNITY_FLASH //flash needs to be able to set these fields +public: +#endif + Vector2f m_Origin; + Vector2f m_Direction; // Direction is always normalized + +public: + Ray2D () {} + Ray2D (const Vector2f& origin, const Vector2f& direction) { m_Origin = origin; SetDirection (direction); } + + const Vector2f& GetOrigin ()const { return m_Origin; } + void SetOrigin (const Vector2f& origin) { m_Origin = origin; } + + const Vector2f& GetDirection () const { return m_Direction; } + void SetDirection (const Vector2f& direction) { AssertIf (!IsNormalized (direction)); m_Direction = direction; } + void SetApproxDirection (const Vector2f& direction) { m_Direction = NormalizeFast (direction); } + + Vector2f GetPoint (const float scale) const { return m_Origin + (scale * m_Direction); } +}; + +#endif diff --git a/Runtime/Geometry/Sphere.cpp b/Runtime/Geometry/Sphere.cpp new file mode 100644 index 0000000..19ddb3c --- /dev/null +++ b/Runtime/Geometry/Sphere.cpp @@ -0,0 +1,57 @@ +#include "UnityPrefix.h" +#include "Sphere.h" + +Sphere MergeSpheres (const Sphere& inSphere0, const Sphere& inSphere1); +float MergeSpheresRadius (const Sphere& inSphere0, const Sphere& inSphere1); + +void Sphere::Set (const Vector3f* inVertices, UInt32 inHowmany) +{ + m_Radius = 0.0F; + m_Center = Vector3f::zero; + UInt32 i; + for (i=0;i<inHowmany;i++) + m_Radius = std::max (m_Radius, SqrMagnitude (inVertices[i])); + m_Radius = sqrt (m_Radius); +} + + +Sphere MergeSpheres (const Sphere& inSphere0, const Sphere& inSphere1) +{ + Vector3f kCDiff = inSphere1.GetCenter() - inSphere0.GetCenter(); + float fLSqr = SqrMagnitude (kCDiff); + float fRDiff = inSphere1.GetRadius() - inSphere0.GetRadius(); + + if (fRDiff*fRDiff >= fLSqr) + return fRDiff >= 0.0 ? inSphere1 : inSphere0; + + float fLength = sqrt (fLSqr); + const float fTolerance = 1.0e-06f; + Sphere kSphere; + + if (fLength > fTolerance) + { + float fCoeff = (fLength + fRDiff) / (2.0 * fLength); + kSphere.GetCenter () = inSphere0.GetCenter () + fCoeff * kCDiff; + } + else + { + kSphere.GetCenter () = inSphere0.GetCenter (); + } + + kSphere.GetRadius () = 0.5F * (fLength + inSphere0.GetRadius () + inSphere1.GetRadius ()); + + return kSphere; +} + +float MergeSpheresRadius (const Sphere& inSphere0, const Sphere& inSphere1) +{ + Vector3f kCDiff = inSphere1.GetCenter() - inSphere0.GetCenter(); + float fLSqr = SqrMagnitude (kCDiff); + float fRDiff = inSphere1.GetRadius () - inSphere0.GetRadius(); + + if (fRDiff*fRDiff >= fLSqr) + return fRDiff >= 0.0 ? inSphere1.GetRadius () : inSphere0.GetRadius (); + + float fLength = sqrt (fLSqr); + return 0.5F * (fLength + inSphere0.GetRadius () + inSphere1.GetRadius ()); +} diff --git a/Runtime/Geometry/Sphere.h b/Runtime/Geometry/Sphere.h new file mode 100644 index 0000000..e19e2f5 --- /dev/null +++ b/Runtime/Geometry/Sphere.h @@ -0,0 +1,79 @@ +#ifndef SPHERE_H +#define SPHERE_H + +#include "Runtime/Serialize/SerializeUtility.h" +#include "Runtime/Math/Vector3.h" +#include <algorithm> +#include "Runtime/Modules/ExportModules.h" + +class Sphere +{ + Vector3f m_Center; + float m_Radius; + + public: + + DECLARE_SERIALIZE (Sphere) + + Sphere () {} + Sphere (const Vector3f& p0, float r) {Set (p0, r);} + + void Set (const Vector3f& p0) {m_Center = p0; m_Radius = 0;} + void Set (const Vector3f& p0, float r) {m_Center = p0; m_Radius = r;} + + void Set (const Vector3f& p0, const Vector3f& p1); + + void Set (const Vector3f* inVertices, UInt32 inHowmany); + + Vector3f& GetCenter () {return m_Center;} + const Vector3f& GetCenter ()const {return m_Center;} + + float& GetRadius () {return m_Radius;} + const float& GetRadius ()const {return m_Radius;} + + bool IsInside (const Sphere& inSphere)const; +}; + +float EXPORT_COREMODULE CalculateSqrDistance (const Vector3f& p, const Sphere& s); +bool Intersect (const Sphere& inSphere0, const Sphere& inSphere1); + + +inline void Sphere::Set (const Vector3f& p0, const Vector3f& p1) +{ + Vector3f dhalf = (p1 - p0) * 0.5; + + m_Center = dhalf + p0; + m_Radius = Magnitude (dhalf); +} + +inline bool Sphere::IsInside (const Sphere& inSphere)const +{ + float sqrDist = SqrMagnitude (GetCenter () - inSphere.GetCenter ()); + if (Sqr (GetRadius ()) > sqrDist + Sqr (inSphere.GetRadius ())) + return true; + else + return false; +} + +inline bool Intersect (const Sphere& inSphere0, const Sphere& inSphere1) +{ + float sqrDist = SqrMagnitude (inSphere0.GetCenter () - inSphere1.GetCenter ()); + if (Sqr (inSphere0.GetRadius () + inSphere1.GetRadius ()) > sqrDist) + return true; + else + return false; +} + +inline float CalculateSqrDistance (const Vector3f& p, const Sphere& s) +{ + return std::max (0.0F, SqrMagnitude (p - s.GetCenter ()) - Sqr (s.GetRadius ())); +} + +template<class TransferFunction> inline +void Sphere::Transfer (TransferFunction& transfer) +{ + TRANSFER (m_Center); + TRANSFER (m_Radius); +} + +#endif diff --git a/Runtime/Geometry/SpriteMeshGenerator.cpp b/Runtime/Geometry/SpriteMeshGenerator.cpp new file mode 100644 index 0000000..ad49bf5 --- /dev/null +++ b/Runtime/Geometry/SpriteMeshGenerator.cpp @@ -0,0 +1,973 @@ +#include "UnityPrefix.h" +#include "SpriteMeshGenerator.h" + +#if ENABLE_SPRITES +#include "Runtime/Profiler/Profiler.h" +#include "Runtime/Graphics/SpriteUtility.h" +#include "Runtime/Math/Vector2.h" +#include "Runtime/Math/FloatConversion.h" +#include "Runtime/Math/Polynomials.h" + +#include "External/libtess2/libtess2/tesselator.h" +#include <queue> + +PROFILER_INFORMATION (gProfileDecompose, "SpriteMeshGenerator.Decompose", kProfilerRender); +PROFILER_INFORMATION (gProfileTraceShape, "SpriteMeshGenerator.TraceShape", kProfilerRender); +PROFILER_INFORMATION (gProfileSimplify, "SpriteMeshGenerator.Simplify", kProfilerRender); + +static const float kHoleAreaLimit = 0.25f; +static const float kMaxTriangles = 1000.0f; +static const float kMaxOverdraw = 4.0f; +static const float kResolution = 960*640*kMaxOverdraw; + +struct edge { + float m_a; + float m_b; + float m_c; + bool m_apos; + bool m_bpos; + edge(){}; + edge(Vector2f p0, Vector2f p1) { + m_a = (p0.y - p1.y); + m_b = (p1.x - p0.x); + m_c = -p0.x*m_a - p0.y*m_b; + + bool aez = (m_a == 0); + bool bez = (m_b == 0); + bool agz = (m_a < 0); + bool bgz = (m_b < 0); + m_apos = aez ? bgz : agz; + m_bpos = bez ? agz : bgz; + } + + float grad(Vector2f p) { return (m_a*p.x + m_b*p.y + m_c); } + int test(Vector2f p) { + float g = grad(p); + return ((g > 0) || ((g == 0) && m_apos)) ? 1 : -1; + } +}; + +static inline int mod(int a, int n) +{ + return a>=n ? a%n : a>=0 ? a : n-1-(-1-a)%n; +} + + +inline float det(const Vector2f& a, const Vector2f& b, const Vector2f& c) +{ + float bax = b.x - a.x; + float acx = a.x - c.x; + float aby = a.y - b.y; + float cay = c.y - a.y; + return (bax * cay) - ( acx * aby ); +} + +inline Vector2f ortho(const Vector2f& v) +{ + return Vector2f(v.y, -v.x); +} + +inline float distance_point_line(Vector2f pq, Vector2f p0, Vector2f p1) +{ + Vector2f v = p1 - p0; + Vector2f w = pq - p0; + + float a = Dot(w, v); + if (a <= 0) + return Magnitude(p0-pq); + + float b = Dot(v, v); + if (b <= a) + return Magnitude(p1-pq); + + float c = a/b; + Vector2f p = p0 + v*c; + return Magnitude(p-pq); +} + +void SpriteMeshGenerator::Decompose(std::vector<Vector2f>* vertices, std::vector<int>* indices) +{ + if (m_paths.size() == 0) + return; + + vertices->clear(); + indices->clear(); + + const int kVertexSize = 2; + const int kPolygonVertices = 3; + + PROFILER_BEGIN(gProfileDecompose, NULL); + TESStesselator* tess = tessNewTess(NULL); + for (std::vector<path>::const_iterator it = m_paths.begin(); it != m_paths.end(); ++it) + { + const path& p = *it; + const std::vector<vertex>& vertices = p.m_path; + if (vertices.size() == 0) + continue; + + tessAddContour(tess, kVertexSize, &vertices[0].p, sizeof(vertex), vertices.size()); + } + int tessError = tessTesselate(tess, TESS_WINDING_NONZERO, TESS_POLYGONS, kPolygonVertices, kVertexSize, NULL); + AssertBreak(tessError == 1); + + const int elemCount = tessGetElementCount(tess); + const TESSindex* elements = tessGetElements(tess); + const TESSreal* real = tessGetVertices(tess); + + for (int e = 0; e < elemCount; ++e) + { + const int* idx = &elements[e * kPolygonVertices]; + + // Extract vertices + for (int i = 0; i < kPolygonVertices; ++i) + { + Assert(idx[i] != TESS_UNDEF); + + float x = real[idx[i]*kVertexSize]; + float y = real[idx[i]*kVertexSize + 1]; + +#define SNAP_VERTEX_POSITION 1 +#if SNAP_VERTEX_POSITION + x = floor(x + 0.5f); + y = floor(y + 0.5f); +#endif +#undef SNAP_VERTEX_POSITION + + Vector2f newVertex(x, y); + + // Reuse vertex + bool reused = false; + for (int v = 0; v < vertices->size(); ++v) + { + const Vector2f& vertex = (*vertices)[v]; + if ((std::abs(vertex.x - x) <= Vector2f::epsilon) && (std::abs(vertex.y - y) <= Vector2f::epsilon)) + { + indices->push_back(v); + reused = true; + break; + } + } + + // New vertex + if (!reused) + { + indices->push_back(vertices->size()); + vertices->push_back(newVertex); + } + } + + // Push polygon + } +#ifdef __MESHGEN_STATS + { + int n=(int)m_paths.size(); + + double mesh_area = 0.0; + double rect_area = 0.0; + for(int i=0; i<n; i++) { + path *p = &m_paths[i]; + int m = p->m_path.size(); + for(int j=0; j<m; j++) { + Vector2f p0 = p->m_path[j].p; + Vector2f p1 = p->m_path[mod(j+1, m)].p; + mesh_area += (p0.x*p1.y)-(p0.y*p1.x); + } + } + mesh_area *= 0.5; + rect_area = m_mask_org.w*m_mask_org.h; + + LogString(Format("Sprite mesh triangle count : %d", ((indices!=NULL) ? (int)(indices->size()/3) : 0))); + LogString(Format("Sprite mesh area (image space) : %.0f", mesh_area)); + LogString(Format("Sprite rect area (image space) : %.0f", rect_area)); + LogString(Format("Sprite diff area (image space) : %.0f", rect_area-mesh_area)); + + } +#endif + + tessDeleteTess(tess); + PROFILER_END +} + +float SpriteMeshGenerator::evaluateLOD(const float areaHint, float area) +{ + // do rough estimation of simplification lod + int triangleCount=0; + int n = (int)m_paths.size(); + + // evaluate optimal triangle count + for(int i=0; i<n; i++) { + path *p = &m_paths[i]; + if (p->isHole()) + triangleCount += 2; + else + triangleCount += (int)p->m_path.size() - 2; + } + float maxTriangleCount = area * areaHint; + float lod = 1.0f-(maxTriangleCount / (float)triangleCount); + return clamp(lod, 0.0f, 1.0f); + +} + +void SpriteMeshGenerator::MakeShape(ColorRGBA32* image, int imageWidth, int imageHeight, float hullTolerance, unsigned char alphaTolerance, bool holeDetection, unsigned int extrude, float bias, int mode) +{ + PROFILER_BEGIN(gProfileTraceShape, NULL); + m_mask_org = mask(image, imageWidth, imageHeight, alphaTolerance, extrude); + m_mask_cur = mask(image, imageWidth, imageHeight, alphaTolerance, extrude); + + std::vector<vertex> outline; + + int sign; + float area; + float areaRect = imageWidth*imageHeight; + float areaTotal = 0; + while(contour(outline, sign, area)) { + if (!holeDetection && sign=='-') + continue; + if (area<(areaRect*kHoleAreaLimit) && sign=='-' && hullTolerance < 0.0f) + continue; + areaTotal += (sign=='+') ? area : -area; + m_paths.push_back(path(outline, imageWidth, imageHeight, sign, area, bias)); + } + PROFILER_END + + PROFILER_BEGIN(gProfileSimplify, NULL); + if (hullTolerance < 0.0f) { + hullTolerance = evaluateLOD(kMaxTriangles/kResolution, areaTotal); + } + // Simplify + for (std::vector<path>::iterator it = m_paths.begin(); it != m_paths.end(); ++it) + { + path& p = *it; + p.simplify(hullTolerance, mode); + } + // Snap + for (std::vector<path>::iterator it = m_paths.begin(); it != m_paths.end(); ++it) + { + path& p = *it; + for (std::vector<vertex>::iterator vit = p.m_path.begin(); vit != p.m_path.end(); ++vit) + { + vertex& vert = *vit; + vert.p.x = Roundf(vert.p.x); + vert.p.y = Roundf(vert.p.y); + } + } + PROFILER_END +} + +inline bool predMinX(const SpriteMeshGenerator::path& a, const SpriteMeshGenerator::path& b) { return a.GetMin().x < b.GetMin().x; } +inline bool predMinY(const SpriteMeshGenerator::path& a, const SpriteMeshGenerator::path& b) { return a.GetMin().y < b.GetMin().y; } +inline bool predMaxX(const SpriteMeshGenerator::path& a, const SpriteMeshGenerator::path& b) { return a.GetMax().x < b.GetMax().x; } +inline bool predMaxY(const SpriteMeshGenerator::path& a, const SpriteMeshGenerator::path& b) { return a.GetMax().y < b.GetMax().y; } + +bool SpriteMeshGenerator::FindBounds(Rectf& bounds) +{ + if (m_paths.size() == 0) + return false; + + const SpriteMeshGenerator::path& minX = *std::min_element(m_paths.begin(), m_paths.end(), predMinX); + const SpriteMeshGenerator::path& minY = *std::min_element(m_paths.begin(), m_paths.end(), predMinY); + const SpriteMeshGenerator::path& maxX = *std::max_element(m_paths.begin(), m_paths.end(), predMaxX); + const SpriteMeshGenerator::path& maxY = *std::max_element(m_paths.begin(), m_paths.end(), predMaxY); + + bounds.x = minX.GetMin().x; + bounds.SetRight(maxX.GetMax().x); + + bounds.y = minY.GetMin().y; + bounds.SetBottom(maxY.GetMax().y); + + return true; +} + +void SpriteMeshGenerator::path::bbox() +{ + float minx=(std::numeric_limits<float>::max)(); + float miny=(std::numeric_limits<float>::max)(); + float maxx=(std::numeric_limits<float>::min)(); + float maxy=(std::numeric_limits<float>::min)(); + + int n=(int)m_path.size(); + for(int i=0; i<n; i++) { + Vector2f p=m_path[i].p; + if (p.x < minx) minx = p.x; + if (p.y < miny) miny = p.y; + if (p.x > maxx) maxx = p.x; + if (p.y > maxy) maxy = p.y; + } + + // clamp to bounds + minx = (minx < 0.0) ? 0.0 : (minx > m_bx) ? m_bx : minx; + miny = (miny < 0.0) ? 0.0 : (miny > m_by) ? m_by : miny; + maxx = (maxx < 0.0) ? 0.0 : (maxx > m_bx) ? m_bx : maxx; + maxy = (maxy < 0.0) ? 0.0 : (maxy > m_by) ? m_by : maxy; + + m_min = Vector2f(minx, miny); + m_max = Vector2f(maxx, maxy); +} + +bool SpriteMeshGenerator::path::dec(int i) +{ + int n = (int)m_path.size(); + if (n<3) + return false; + + Vector2f a = m_path[mod(i-1, n)].p; + Vector2f b = m_path[mod(i+0, n)].p; + Vector2f c = m_path[mod(i+1, n)].p; + Vector2f ab = a-b; + Vector2f bc = b-c; + Vector2f na = NormalizeSafe(Vector2f(-ab.y, ab.x)); + Vector2f nb = NormalizeSafe(Vector2f(-bc.y, bc.x)); + Vector2f no = NormalizeSafe(nb+na); + m_path[mod(i, n)].n = no; + return true; +} + +bool SpriteMeshGenerator::path::inf(int i) +{ + int n = (int)m_path.size(); + if (n<3) + return false; + Vector2f a = m_path[mod(i-1, n)].p; + Vector2f b = m_path[mod(i+0, n)].p; + Vector2f c = m_path[mod(i+1, n)].p; + m_path[i].s = edge(a,c).test(b); + return true; +} + +static int dir(Vector2f p0, Vector2f p1) +{ + int di[3*3] = { + 0, 1, 2, + 7, -1, 3, + 6, 5, 4 + }; + Vector2f dt = p0 - p1; + int dx = (dt.x > 0.0f) ? 1 : (dt.x < 0.0f) ? -1 : 0; + int dy = (dt.y > 0.0f) ? 1 : (dt.y < 0.0f) ? -1 : 0; + int d = 4 + 3*dx - dy; + return (d>=0 || d<=8) ? di[d] : -1; +} + +static bool min_positive(float a, float b, float& res) +{ + if(a > 0 && b > 0) + res = a < b ? a : b; + else + res = a > b ? a : b; + return ((res > 0) || CompareFloatRobust(res, 0.0)); +} + +#define LE 0x1 +#define RE 0x2 +#define BE 0x4 +#define TE 0x8 + +bool SpriteMeshGenerator::path::clip_test(Vector2f p, int side) +{ + switch( side ) { + case LE: return p.x >= m_min.x; + case RE: return p.x <= m_max.x; + case TE: return p.y >= m_min.y; + case BE: return p.y <= m_max.y; + } + return false; +} + +Vector2f SpriteMeshGenerator::path::clip_isec(Vector2f p, Vector2f q, int e) +{ + double a = (q.y - p.y) / (q.x - p.x); + double b = p.y - p.x * a; + double x, y; + switch(e) { + case LE: + case RE: + x = (e == LE) ? m_min.x : m_max.x; + y = x * a + b; + break; + case TE: + case BE: + y = (e == TE) ? m_min.y : m_max.y; + x = (IsFinite(a)) ? (y - b) / a : p.x; + break; + } + return Vector2f(x,y); +} + +void SpriteMeshGenerator::path::clip_edge(int e) +{ + int n=(int)m_path.size(); + std::vector<vertex> cpath; + + for (int i=0 ; i<n; i++) { + Vector2f s = m_path[mod(i+0, n)].p; + Vector2f p = m_path[mod(i+1, n)].p; + Vector2f c; + if (clip_test(p, e)) { + if (!clip_test(s, e) ) { + c = clip_isec(p, s, e); + cpath.push_back(vertex(c)); + } + cpath.push_back(vertex(p)); + } + else + if (clip_test(s, e)) { + c = clip_isec(s, p, e); + cpath.push_back(vertex(c)); + } + } + m_path.clear(); + m_path=cpath; +} + +void SpriteMeshGenerator::path::clip() +{ + clip_edge(LE); + clip_edge(RE); + clip_edge(TE); + clip_edge(BE); +} + +static bool lseg_intersect(Vector2f p1, Vector2f p2, Vector2f p3, Vector2f p4) +{ + Vector2f e43 = p4-p3; + Vector2f e21 = p2-p1; + Vector2f e13 = p1-p3; + + double dn = e43.y*e21.x - e43.x*e21.y; + double na = e43.x*e13.y - e43.y*e13.x; + double nb = e21.x*e13.y - e21.y*e13.x; + + // coincident? + if ((fabs(na) < Vector2f::epsilon) && + (fabs(nb) < Vector2f::epsilon) && + (fabs(dn) < Vector2f::epsilon)) { + return false; + } + + // parallel ? + if (fabs(dn) < Vector2f::epsilon) + return false; + + // collinear ? + double mua = na / dn; + double mub = nb / dn; + if ((mua < 0) || (mua > 1) || (mub < 0) || (mub > 1)) + return false; + + return true; +} + +int SpriteMeshGenerator::path::self_intersect(Vector2f p0, Vector2f p1) +{ + int n=(int)m_path.size(); + for(int i=0; i<n; i++) { + Vector2f p2 = m_path[mod(i+0, n)].p; + Vector2f p3 = m_path[mod(i+1, n)].p; + if ((p2==p0) || + (p3==p1) || + (p1==p2) || + (p0==p3)) + continue; + if (lseg_intersect(p0, p1, p2, p3)) + return 1; + } + return 0; +} + +bool SpriteMeshGenerator::path::cvx_cost(int i) +{ + int n = (int)m_path.size(); + if (n<5) + return false; + + vertex *v = &m_path[i]; + Vector2f sn = m_path[mod(i-1,n)].n; + Vector2f tn = m_path[mod(i+1,n)].n; + // detect high convexity -> better triangulation for cyclic geometric shape + float q = Dot(sn,tn); + if ((q < 0.000) || + CompareFloatRobust(q, 0.0) || + CompareFloatRobust(q, 1.0) ) { + v->cost=s_cost(-1, 0.0); + return true; + } + + Vector2f s0 = m_path[mod(i-1,n)].p; + Vector2f t0 = m_path[mod(i+1,n)].p; + Vector2f p0 = v->p; + Vector2f a0 = ortho(sn); + Vector2f c0 = ortho(p0-s0); + Vector2f c1 = s0-t0; + Vector2f b0 = tn-sn; + float c = Dot(c0, c1); + float a = Dot(tn, a0); + float b = Dot(b0, c0) + Dot(c1, a0); + + float x0=-1,x1=-1,w,d; + if (QuadraticPolynomialRootsGeneric(a, -b, c, x0, x1) && min_positive(x0, x1, w)) { + Vector2f r0 = m_path[mod(i-2, n)].p; + Vector2f u0 = m_path[mod(i+2, n)].p; + Vector2f s1 = s0 + sn*w; + Vector2f t1 = t0 + tn*w; + + float d0 = det(r0, s1, s0); + float d1 = det(s1, p0, s0); + float d2 = det(p0, t1, t0); + float d3 = det(t1, u0, t0); + float d = d0+d1+d2+d3; + v->cost = s_cost(d, w); + return true; + } + else { + // this should not ever happen.. but handle it anyway + d = -1.0; + w = 0.0; + v->cost = s_cost(d, w); + return false; + } +} + +bool SpriteMeshGenerator::path::cve_cost(int i) +{ + int n=(int)m_path.size(); + if (n<3) + return 0; + vertex *v = &m_path[i]; + Vector2f s = m_path[mod(i-1,n)].p; + Vector2f t = m_path[mod(i+1,n)].p; + Vector2f u = v->p; + float d = det(s,t,u); + if ((d>0) || CompareFloatRobust(d, 0.0)) { + v->cost = s_cost(d, 0.0); + return true; + } + else { + v->cost = s_cost(-1.0, 0.0); + return false; + } +} + +int SpriteMeshGenerator::path::min_cost() +{ + int n = (int)m_path.size(); + int min_i = -1; + float min_c = (std::numeric_limits<float>::max)(); + for (int i=0; i<n; i++) { + vertex v = m_path[i]; + if (v.cost.c < 0) + continue; + float c = v.cost.c + v.c; + if (c < min_c) { + min_c = c; + min_i = i; + } + } + return min_i; +} + +bool SpriteMeshGenerator::path::select() +{ + int n = (int)m_path.size(); + if (n<5) + return false; + + int i = 0; + int m = 0; + bool found = false; + do { + i = min_cost(); + if (i<0) + return false; + struct vertex *v = &m_path[i]; + struct vertex *s = &m_path[mod(i-1, n)]; + struct vertex *t = &m_path[mod(i+1, n)]; + struct s_cost cost = v->cost; + int isec = 0; + if (v->s > 0) { + Vector2f s0 = s->p; + Vector2f t0 = t->p; + // check self intersection + isec = self_intersect(s0, t0); + + if (isec) v->cost.c = -1; + else { + s->c += cost.c; + t->c += cost.c; + m_invalid.push_back(mod(i+0, n)); + m_invalid.push_back(mod(i+1, n)); + found = true; + } + } + else { + struct vertex *r = &m_path[mod(i-2, n)]; + struct vertex *u = &m_path[mod(i+2, n)]; + Vector2f r0 = r->p; + Vector2f u0 = u->p; + Vector2f s0 = s->p + s->n*cost.w; + Vector2f t0 = t->p + t->n*cost.w; + // check self intersection + isec |= self_intersect(r0, s0); + isec |= self_intersect(s0, t0); + isec |= self_intersect(t0, u0); + + if (isec) v->cost.c = -1; + else { + s->p = s0; + t->p = t0; + s->c += cost.c; + t->c += cost.c; + m_invalid.push_back(mod(i-2, n)); + m_invalid.push_back(mod(i-1, n)); + m_invalid.push_back(mod(i+1, n)); + m_invalid.push_back(mod(i+2, n)); + found = true; + } + } + }while((m++ < n) && !found); + + if (found) { + m_path.erase(m_path.begin()+i); + // fix invalid indices after erase + m = (int)m_invalid.size(); + for (int k=0; k<m; k++) { + if (m_invalid[k] > i) + m_invalid[k] = m_invalid[k]-1; + } + } + return found; +} + +int SpriteMeshGenerator::path::find_max_distance(int i0) +{ + int n = m_path.size(); + + Vector2f a = m_path[i0].p; + float dm = -1; + int mi = -1; + + for (int i=0; i<n; i++) { + Vector2f b = m_path[mod(i, n)].p; + float ba = Magnitude(b-a); + if (ba < dm) + continue; + dm = ba; + mi = i; + } + return mi; +} + +int SpriteMeshGenerator::path_segment::max_distance(std::vector<vertex> path, int i0, int i1) +{ + int n = path.size(); + + Vector2f a = path[i0].p; + Vector2f b = path[i1].p; + + float dm = -1; + int mi = -1; + m_cnt=0; + for (int i=i0; i != i1; i=mod(++i, n), m_cnt++) { + float dq = distance_point_line(path[i].p, a, b); + if (dq < dm) + continue; + dm = dq; + mi = i; + } + return mi; +} + +void SpriteMeshGenerator::path::simplify(float q, int mode) +{ + m_path.clear(); + m_path = m_path0; + int m; + int n = (int)m_path.size(); + int lim = (float)n*(1.0f - clamp(q, 0.0f, 1.0f)); + + if (n < 5) goto bail_out; + + if (mode==kPathEmbed) { + if (lim < 5) lim=5; + + // mark all vertices invalid + for (int i=0; i<n; i++) + m_invalid.push_back(i); + + do { + n = (int)m_path .size(); + m = (int)m_invalid.size(); + for (int i=0; i<m; i++) { + dec(m_invalid[i]); + inf(m_invalid[i]); + } + for (int i=0; i<m; i++) { + int k = m_invalid[i]; + if (m_path[k].s > 0) + cve_cost(k); + else + cvx_cost(k); + } + m_invalid.clear(); + if (select() == false) + break; + }while(n > lim); + } + else { + if (lim < 4) lim=4; + + int i0 = find_max_distance( 0 ); + int i1 = find_max_distance(i0 ); + + path_segment ls = path_segment(m_path, i0, i1); + path_segment rs = path_segment(m_path, i1, i0); + + std::priority_queue<path_segment, std::vector<path_segment>, compare_path_segment> pq; + + if (ls.m_mx>=0) pq.push(ls); + if (rs.m_mx>=0) pq.push(rs); + + std::vector<bool> select(n); + std::fill(select.begin(), select.end(), false); + + select[i0] = true; + select[i1] = true; + + int count=2; + while (!pq.empty()) { + path_segment ts = pq.top(); + pq.pop(); + + select[ts.m_mx]=true; + if (++count == lim) + break; + // split + ls = path_segment(m_path, ts.m_i0, ts.m_mx); + if (ls.m_mx >= 0) pq.push(ls); + rs = path_segment(m_path, ts.m_mx, ts.m_i1); + if (rs.m_mx >= 0) pq.push(rs); + } + + m_path.clear(); + for (int i=0; i<n; i++) { + if (select[i]==1) + m_path.push_back(m_path0[i]); + } + } +bail_out: + if (m_sign == '+' && mode==1) + clip(); +} + +void SpriteMeshGenerator::path::fit(std::vector<int>& ci, int i0, int i1) +{ + int n = (int)m_path.size(); + + if ((mod(i0+1, n) == i1) || (i0==i1)) { + ci.push_back(i1); + return; + } + + Vector2f a = m_path[i0].p; + Vector2f b = m_path[i1].p; + edge e = edge(a, b); + int im = -1; + float qm = -1; + int ic = i0; + do { + float qc = fabs(e.grad(m_path[ic].p)); + if (qc > qm) { + im = ic; + qm = qc; + } + if (ic == i1) + break; + ic = mod(ic+1, n); + }while(1); + + float lim = std::max(fabs(e.m_a)*0.5, fabs(e.m_b)*0.5); + if ( (qm <= lim) || (im < 0)) { + ci.push_back(i1); + return; + } + fit(ci, i0, im); + fit(ci, im, i1); +} + +bool SpriteMeshGenerator::path::opt(float bias) +{ + int n = (int)m_path.size(); + if (n<3) + return false; + + std::vector<int> cp; + std::vector<int> ci; + + int s = -1; + int dt[8] = {0}; + cp.push_back(0); + for (int i=0; i<n; i++) { + Vector2f p0 = m_path[i].p; + Vector2f p1 = m_path[mod(i+1, n)].p; + + int d = dir(p0, p1); + if (d < 0) + continue; + dt[d] = 1; + if (s < 0) { + s = d; + continue; + } + // cut path, if direction change is not possible for straight line + if (!(d == mod(s-1, 8) || d == mod(s+1, 8) || d == s) || + ((dt[0] + dt[1] + + dt[2] + dt[3] + + dt[4] + dt[5] + + dt[6] + dt[7] ) > 2)) { + memset(dt, 0, 8*sizeof(int)); + s = -1; + cp.push_back(i); + } + } + // fit sub paths to straight line + int m = (int)cp.size(); + for (int i=0; i<m; i++) + fit(ci, cp[i], cp[mod(i+1, m)]); + + //rm extra vertices + std::vector<vertex> tmp = m_path; + m_path.clear(); + for(int i=0; i<ci.size(); i++) + m_path.push_back(tmp[ci[i]]); + + // unit normals + n = (int)m_path.size(); + for (int i=0; i<n; i++) + dec(i); + + // bias + for (int i=0; i<n; i++) + m_path[i].p += m_path[i].n*bias; + + return true; +} + +bool SpriteMeshGenerator::invmask(std::vector<vertex>& outline) +{ + int n = (int)outline.size(); + if (n <= 0) + return false; + + int xa = (int)outline[ 0].p.x; + Vector2f pp = outline[n-1].p; + + for (int i=0; i<n; i++) { + Vector2f p0 = outline[i].p; + + while (((i+1)<n) && (p0.y == outline[i+1].p.y) ) { + int d = dir(p0, outline[i+1].p); + if ((d==1 && (pp.y < p0.y)) || + (d==5 && (pp.y > p0.y)) ) + p0 = outline[i+1].p; + i++; + } + int y = (int)p0.y; + int x0 = min(xa, (int)p0.x); + int x1 = max(xa, (int)p0.x); + for (int x=x0; x<x1; x++) + m_mask_cur.inv(x, y); + + if (((i+1) < n) && + (pp.y != p0.y) && (outline[i+1].p.y == pp.y) ) { + y = (int)p0.y; + x0 = min(xa, (int)p0.x); + x1 = max(xa, (int)p0.x); + for (int x=x0; x<x1; x++) + m_mask_cur.inv(x, y); + } + pp = p0; + } + for (int i=0; i<n; i++) { + Vector2f p = outline[i].p; + m_mask_cur.rst(p.x, p.y); + } + return true; +} + +bool SpriteMeshGenerator::trace(Vector2f p0, Vector2f p1, Vector2f &p) +{ + static int dt[8][2] = { + { -1, 0 }, + { -1, -1 }, + { 0, -1 }, + { 1, -1 }, + { 1, 0 }, + { 1, 1 }, + { 0, 1 }, + { -1, 1 } + }; + + int t0 = dir(p0, p1); + if (t0 < 0) + goto error; + + for (int i = 0; i < 8; i++) { + int t = (t0 + i) % 8; + int x = (int)p1.x + dt[t][0]; + int y = (int)p1.y + dt[t][1]; + if (m_mask_cur.tst(x, y)) { + p = Vector2f(x, y); + return true; + } + } +error: + p = Vector2f(-1, -1); + return false; +} + +bool SpriteMeshGenerator::contour(std::vector<vertex>& outline, int &sign, float &area) +{ + do { + outline.clear(); + int b = m_mask_cur.first(); + if (b < 0) + return false; + + int x = b % m_mask_cur.w; + int y = b / m_mask_cur.w; + Vector2f curr = Vector2f (x, y); + Vector2f stop; + Vector2f prev; + Vector2f next; + area = 0.0; + sign = m_mask_org.tst(x, y) ? '+' : '-'; + stop = curr; + next = curr; + curr.x = curr.x-1; + + do { + prev = curr; + curr = next; + outline.push_back(vertex(curr)); + if (trace(prev, curr, next) == false) + break; + + area += curr.x * next.y - + curr.y * next.x; + if (next == stop) + break; + } while(true); + + invmask(outline); + if (fabs(area)<4) { + area=0; + continue; + } + if (((sign=='+') && (area < 0)) || + ((sign=='-') && (area > 0)) ) + std::reverse(outline.begin(), outline.end()); + area = fabs(area); + break; + }while(1); + return true; +} +#endif //ENABLE_SPRITES diff --git a/Runtime/Geometry/SpriteMeshGenerator.h b/Runtime/Geometry/SpriteMeshGenerator.h new file mode 100644 index 0000000..927b629 --- /dev/null +++ b/Runtime/Geometry/SpriteMeshGenerator.h @@ -0,0 +1,269 @@ +#pragma once +#include "Configuration/UnityConfigure.h" + +#if ENABLE_SPRITES + +#include "Runtime/BaseClasses/NamedObject.h" +#include "Runtime/Math/Rect.h" +#include "Runtime/Serialize/TransferFunctions/SerializeTransfer.h" +#include "Runtime/Shaders/Material.h" +#include "Runtime/Graphics/Texture2D.h" +#include "Runtime/Utilities/dynamic_bitset.h" + + +class SpriteMeshGenerator +{ +public: + struct s_cost + { + s_cost(){}; + s_cost(float _c, float _w) + { + c=_c; + w=_w; + }; + float c; + float w; + }; + + struct vertex + { + vertex(){}; + vertex(Vector2f pos) + { + p = pos; + }; + Vector2f p; + Vector2f n; // normal + int s; // sign -> indicating concavity + float c; + struct s_cost cost; + }; + + class path_segment + { + public: + path_segment(std::vector<vertex> path, int i0, int i1) + { + m_i0 = i0; + m_i1 = i1; + m_mx = max_distance(path, i0, i1); + }; + + int m_i0; + int m_i1; + int m_mx; + int m_cnt; + private: + int max_distance(std::vector<vertex> path, int i0, int i1); + }; + + class compare_path_segment { + public: + bool operator()(path_segment& s0, path_segment& s1) + { + return (s0.m_cnt < s1.m_cnt); + } + }; + + class path + { + public: + path(){}; + path(const std::vector<vertex>& p, int w, int h, int sign, float area, float bias) + { + m_bx = w; + m_by = h; + m_area = area; + m_sign = sign; + m_path = p; + opt (bias); + bbox(); + m_path0 = m_path; + } + + std::vector<vertex> m_path; + + void bbox (); + void simplify (float q, int mode); + + bool isHole() const { return m_sign == '-'; } + + const Vector2f& GetMin() const { return m_min; } + const Vector2f& GetMax() const { return m_max; } + + private: + int find_max_distance(int i0); + + void fit (std::vector<int>& ci, int i0, int i1); + bool opt (float bias); + + bool dec (int i); + bool inf (int i); + bool select(); + bool cvx_cost(int i); + bool cve_cost(int i); + int min_cost(); + int self_intersect(Vector2f p0, Vector2f p1); + + void clip(); + void clip_edge(int e); + bool clip_test(Vector2f p, int side); + Vector2f clip_isec(Vector2f p, Vector2f q, int e); + + int m_bx; + int m_by; + int m_sign; + float m_area; + Vector2f m_min; + Vector2f m_max; + std::vector<vertex> m_path0; + std::vector<int> m_invalid; + }; + + struct mask + { + mask(){}; + mask(ColorRGBA32 *img, int width, int height, unsigned char acut, unsigned int dn) + { + w = width; + h = height; + int n = w*h; + dynamic_bitset bv; + bv.resize(n); + for (int y=0; y<height; y++) + for (int x=0; x<width; x++) { + if (img[x+y*width].a > acut) + bv.set(x+y*w); + } + + if (dn) + this->dilate(dn, bv); + + w+=1; + h+=1; + m_bv.resize(w*h); + for (int y=0; y<height; y++) + for (int x=0; x<width; x++) { + if (bv.test(x+y*width)) { + m_bv.set((x+0)+(y+0)*w); + m_bv.set((x+1)+(y+1)*w); + m_bv.set((x+0)+(y+1)*w); + m_bv.set((x+1)+(y+0)*w); + } + } + } + + bool tst(int x, int y) + { + if ((x<0) || (x>=w) || + (y<0) || (y>=h) ) + return 0; + int i=x+y*w; + return m_bv.test(i); + } + + void set(int x, int y) + { + if ((x<0) || (x>=w) || + (y<0) || (y>=h) ) + return; + int i=x+y*w; + m_bv.set(i); + } + + void rst(int x, int y) + { + if ((x<0) || (x>=w) || + (y<0) || (y>=h) ) + return; + int i=x+y*w; + m_bv[i]=0; + } + + void inv(int x, int y) + { + if ((x<0) || (x>=w) || + (y<0) || (y>=h) ) + return; + int i=x+y*w; + m_bv[i].flip(); + } + + int first() + { + int n=(int)m_bv.size(); + for (int i=0; i<n; i++) + if (m_bv.test(i)) + return i; + return -1; + } + + bool dilate(int n, dynamic_bitset &bv) + { + if ((w==0) || (h==0)) + return false; + UInt32 *md = new UInt32[w*h]; + if (!mdist(md, bv)) + return false; + + for (int i=0; i<w*h; i++) { + if (md[i] <= n) + bv.set(i); + } + delete md; + return true; + } + int w; + int h; + dynamic_bitset m_bv; + private: + bool mdist(UInt32 *md, dynamic_bitset& bv ) + { + if (!md) + return false; + + for (int y=0; y<h; y++) + for (int x=0; x<w; x++) { + int i = x+y*w; + if (bv.test(i)) + md[i] = 0; + else { + md[i] = w+h; + if (y>0) md[i] = min(md[i], md[i-w]+1); + if (x>0) md[i] = min(md[i], md[i-1]+1); + } + } + + for (int y=h-1; y>=0; y--) + for (int x=w-1; x>=0; x--) { + int i = x+y*w; + if ((y+1) < h) md[i] = min(md[i], md[i+w]+1); + if ((x+1) < w) md[i] = min(md[i], md[i+1]+1); + } + return true; + } + }; + +public: + void Decompose(std::vector<Vector2f>* vertices, std::vector<int>* indices); + void MakeShape(ColorRGBA32* image, int imageWidth, int imageHeight, float hullTolerance, unsigned char alphaTolerance, bool holeDetection, unsigned int extrude, float bias, int mode); + bool FindBounds(Rectf& bounds); + + const std::vector<path>& GetPaths() const { return m_paths; } + + +private: + bool trace(Vector2f p0, Vector2f p1, Vector2f &p); + bool invmask(std::vector<vertex>& outline); + bool contour(std::vector<vertex>& outline, int &sign, float &area); + + std::vector<path> m_paths; + + float evaluateLOD(const float areaHint, float area); + + struct mask m_mask_org; + struct mask m_mask_cur; +}; + +#endif //ENABLE_SPRITES diff --git a/Runtime/Geometry/TangentSpaceCalculation.cpp b/Runtime/Geometry/TangentSpaceCalculation.cpp new file mode 100644 index 0000000..f7b17f1 --- /dev/null +++ b/Runtime/Geometry/TangentSpaceCalculation.cpp @@ -0,0 +1,534 @@ +#include "UnityPrefix.h" +#include "TangentSpaceCalculation.h" +#include "Runtime/Math/Matrix3x3.h" +#include "Runtime/Math/Vector2.h" +#include "Plane.h" +#include <vector> + +using std::vector; +/* +void CreateTangentSpaceTangents(const Vector3f* vertex, const Vector2f* texcoord, const Vector3f* normals, Matrix3x3f* avgOrthonormalBases, int vertexCount, const UInt16* indices, int triangleCount) +{ + Vector3f *tan1 = new Vector3f[vertexCount * 2]; + Vector3f *tan2 = tan1 + vertexCount; + memset(tan1, 0, vertexCount * sizeof(Vector3f) * 2); + + for (int a = 0; a < triangleCount; a++) + { + int i1 = indices[0]; + int i2 = indices[1]; + int i3 = indices[2]; + + const Vector3f& v1 = vertex[i1]; + const Vector3f& v2 = vertex[i2]; + const Vector3f& v3 = vertex[i3]; + + const Vector2f& w1 = texcoord[i1]; + const Vector2f& w2 = texcoord[i2]; + const Vector2f& w3 = texcoord[i3]; + + float x1 = v2.x - v1.x; + float x2 = v3.x - v1.x; + float y1 = v2.y - v1.y; + float y2 = v3.y - v1.y; + float z1 = v2.z - v1.z; + float z2 = v3.z - v1.z; + + float s1 = w2.x - w1.x; + float s2 = w3.x - w1.x; + float t1 = w2.y - w1.y; + float t2 = w3.y - w1.y; + + float r = 1.0F / (s1 * t2 - s2 * t1); + Vector3f sdir((t2 * x1 - t1 * x2) * r, (t2 * y1 - t1 * y2) * r, (t2 * z1 - t1 * z2) * r); + Vector3f tdir((s1 * x2 - s2 * x1) * r, (s1 * y2 - s2 * y1) * r, (s1 * z2 - s2 * z1) * r); + + tan1[i1] += sdir; + tan1[i2] += sdir; + tan1[i3] += sdir; + + tan2[i1] += tdir; + tan2[i2] += tdir; + tan2[i3] += tdir; + + indices += 3; + } + + + for (long a = 0; a < vertexCount; a++) + { + Vector3f normal = normals[a]; + Vector3f tangent = tan1[a]; + + + // Gram-Schmidt orthogonalize + (tangent - normal * (normal * tangent))); + + // Calculate handedness + //tangent[a].w = (normal % tangent * tan2[a] < 0.0F) ? -1.0F : 1.0F; + + avgOrthonormalBases[a].SetOrthoNormalBasis (tangent, Cross (normal, tangent), normal); + } + + delete[] tan1; +}*/ + +/* +void CreateTangentSpaceTangents (const Vector3f* pPositions, const Vector2f* tex, const Vector3f* normals, Matrix3x3f* avgOrthonormalBases, int vertexCount, const UInt16* indices, int faceCount) +{ + vector<Vector3f> sVector, tVector; + vector<Vector3f> avgS, avgT; + + avgS.resize (vertexCount); + avgT.resize (vertexCount); + + sVector.reserve (faceCount * 3); + tVector.reserve (faceCount * 3); + + // for each face, calculate its S,T & SxT vector, & store its edges + for (int f = 0; f < faceCount * 3; f += 3 ) + { + Vector3f edge0; + Vector3f edge1; + Vector3f s; + Vector3f t; + + // create an edge out of x, s and t + edge0.x = pPositions[ indices[ f + 1 ] ].x - pPositions[ indices[ f ] ].x; + edge0.y = tex[ indices[ f + 1 ] ].x - tex[ indices[ f ] ].x; + edge0.z = tex[ indices[ f + 1 ] ].y - tex[ indices[ f ] ].y; + + // create an edge out of x, s and t + edge1.x = pPositions[ indices[ f + 2 ] ].x - pPositions[ indices[ f ] ].x; + edge1.y = tex[ indices[ f + 2 ] ].x - tex[ indices[ f ] ].x; + edge1.z = tex[ indices[ f + 2 ] ].y - tex[ indices[ f ] ].y; + + Vector3f sxt = Cross (edge0, edge1); + + float a = sxt.x; + float b = sxt.y; + float c = sxt.z; + + float ds_dx = 0.0F; + if ( Abs( a ) > Vector3f::epsilon ) + ds_dx = -b / a; + + float dt_dx = 0.0F; + if ( Abs( a ) > Vector3f::epsilon ) + dt_dx = -c / a; + + // create an edge out of y, s and t + edge0.x = pPositions[ indices[ f + 1 ] ].y - pPositions[ indices[ f ] ].y; + // create an edge out of y, s and t + edge1.x = pPositions[ indices[ f + 2 ] ].y - pPositions[ indices[ f ] ].y; + + sxt = Cross (edge0, edge1); + + a = sxt.x; + b = sxt.y; + c = sxt.z; + + float ds_dy = 0.0F; + if ( Abs( a ) > Vector3f::epsilon ) + ds_dy = -b / a; + + float dt_dy = 0.0F; + if ( Abs( a ) > Vector3f::epsilon ) + dt_dy = -c / a; + + // create an edge out of z, s and t + edge0.x = pPositions[ indices[ f + 1 ] ].z - pPositions[ indices[ f ] ].z; + // create an edge out of z, s and t + edge1.x = pPositions[ indices[ f + 2 ] ].z - pPositions[ indices[ f ] ].z; + + sxt = Cross (edge0, edge1); + + a = sxt.x; + b = sxt.y; + c = sxt.z; + + float ds_dz = 0.0F; + if ( Abs( a ) > Vector3f::epsilon ) + ds_dz = -b / a; + + float dt_dz = 0.0F; + if ( Abs( a ) > Vector3f::epsilon ) + dt_dz = -c / a; + + // generate coordinate frame from the gradients + s = Vector3f( ds_dx, ds_dy, ds_dz ); + t = Vector3f( dt_dx, dt_dy, dt_dz ); + + s = NormalizeSafe (s); + t = NormalizeSafe (t); + + // save vectors for this face + sVector.push_back( s ); + tVector.push_back( t ); + } + + + for ( int p = 0; p < vertexCount; p ++ ) + { + avgS[p] = Vector3f::zero; + avgT[p] = Vector3f::zero; + } + + // go through faces and add up the bases for each vertex + for ( int face = 0; face < faceCount; ++face ) + { + // sum bases, so we smooth the tangent space across edges + avgS[ indices[ face * 3 ] ] += sVector[ face ]; + avgT[ indices[ face * 3 ] ] += tVector[ face ]; + + avgS[ indices[ face * 3 + 1 ] ] += sVector[ face ]; + avgT[ indices[ face * 3 + 1 ] ] += tVector[ face ]; + + avgS[ indices[ face * 3 + 2 ] ] += sVector[ face ]; + avgT[ indices[ face * 3 + 2 ] ] += tVector[ face ]; + } + + // now renormalize + for ( int p = 0; p < vertexCount; p ++ ) + { + Vector3f normal = normals[p]; + + + //OrthoNormalize (&normal, &avgS[p], &avgT[p]); + + avgOrthonormalBases[p].SetOrthoNormalBasis (avgS[p], avgT[p], normal); + + #if DEBUGMODE + float det = avgOrthonormalBases[p].GetDeterminant (); + AssertIf (!CompareApproximately (det, 1.0F,0.001) && !CompareApproximately (det, -1.0F,0.001)); + #endif +// AssertIf (!CompareApproximately (avgOrthonormalBases[p].MultiplyPoint3Transpose (Vector3f (0,0,1)), normal)); + } +}*/ + +/* +void CreateTangentSpaceTangents (const Vector3f* pPositions, const Vector2f* tex, const Vector3f* normals, Matrix3x3f* avgOrthonormalBases, int vertexCount, const UInt16* indices, int faceCount) +{ + + + Vector3f v0,v1,v2, tanu, tanv; + + Vector3f p0,p1,p2; + + Vector3f d1,d2; + + float uv[3][2],det,u,v,l1,l2; + + int i,j,k; + vector<Vector3f> sVector, tVector; + vector<Vector3f> avgS, avgT; + + avgS.resize (vertexCount); + avgT.resize (vertexCount); + + sVector.reserve (faceCount * 3); + tVector.reserve (faceCount * 3); + + for( i=0;i<vertexCount;i++ ) + + { + + sVector[i]=Vector3f(0.0, 0.0, 0.0); + + tVector[i]=Vector3f(0.0, 0.0, 0.0); + + } + + k=0; + for( i=0;i<faceCount;i++,k+=3 ) + + { + + v0=*((Vector3f *)&pPositions[indices[k]]); + + v1=*((Vector3f *)&pPositions[indices[k+1]])-v0; + + v2=*((Vector3f *)&pPositions[indices[k+2]])-v0; + + uv[0][0]=-tex[indices[k]].x; + + uv[0][1]=tex[indices[k]].y; + + uv[1][0]=-tex[indices[k+1]].x-uv[0][0]; + + uv[1][1]=tex[indices[k+1]].y-uv[0][1]; + + uv[2][0]=-tex[indices[k+2]].x-uv[0][0]; + + uv[2][1]=tex[indices[k+2]].y-uv[0][1]; + + det=(uv[1][0]*uv[2][1])-(uv[2][0]*uv[1][1]); + + + + if (fabsf(det)<0.000000001f){ + continue; + } + + u=0; v=0; + + u-=uv[0][0]; v-=uv[0][1]; + + p0=v0+v1*((u*uv[2][1]-uv[2][0]*v)/det)+v2*((uv[1][0]*v-u*uv[1][1])/det); + + + + u=1; v=0; + + u-=uv[0][0]; v-=uv[0][1]; + + p1=v0+v1*((u*uv[2][1]-uv[2][0]*v)/det)+v2*((uv[1][0]*v-u*uv[1][1])/det); + + + + u=0; v=1; + + u-=uv[0][0]; v-=uv[0][1]; + + p2=v0+v1*((u*uv[2][1]-uv[2][0]*v)/det)+v2*((uv[1][0]*v-u*uv[1][1])/det); + + + + d1=p2-p0; + + d2=p1-p0; + + l1=Magnitude(d1); + + l2=Magnitude(d2); + + d1*=1.0f/l1; + + d2*=1.0f/l2; + + + + j=indices[k]; + + sVector[j].x+=d1.x; sVector[j].y+=d1.y; sVector[j].z+=d1.z; + + tVector[j].x+=d2.x; tVector[j].y+=d2.y; tVector[j].z+=d2.z; + + + + j=indices[k+1]; + + sVector[j].x+=d1.x; sVector[j].y+=d1.y; sVector[j].z+=d1.z; + + tVector[j].x+=d2.x; tVector[j].y+=d2.y; tVector[j].z+=d2.z; + + + + j=indices[k+2]; + + sVector[j].x+=d1.x; sVector[j].y+=d1.y; sVector[j].z+=d1.z; + + tVector[j].x+=d2.x; tVector[j].y+=d2.y; tVector[j].z+=d2.z; + + } + + + + for( i=0;i<vertexCount;i++ ) + + { + + //v0.vec(vert[i].tanu[0],vert[i].tanu[1],vert[i].tanu[2]); + v0 = Vector3f(sVector[i].x,sVector[i].y,sVector[i].z); + //v0.normalize(); + v0 = NormalizeRobust(v0); + + //v1.vec(vert[i].tanv[0],vert[i].tanv[1],vert[i].tanv[2]); + v1 = Vector3f(tVector[i].x,tVector[i].y,tVector[i].z); + //v1 = NormalizeRobust(v1); + + Vector3f n(normals[i].x,normals[i].y,normals[i].z); + + if (SqrMagnitude(v1)<0.0001f) + { + v1 = Cross(n,v0); + } + + v1 = NormalizeRobust(v1); + + + sVector[i].x=v0.x; + + sVector[i].y=v0.y; + + sVector[i].z=v0.z; + + tVector[i].x=v1.x; + + tVector[i].y=v1.y; + + tVector[i].z=v1.z; + + avgOrthonormalBases[i].SetOrthoNormalBasis (tVector[i], sVector[i], n); + } + +}*/ + +void CreateTangentSpaceTangents (const Vector3f* pPositions, const Vector2f* tex, const Vector3f* normals, Matrix3x3f* avgOrthonormalBases, int vertexCount, const int* indices, int faceCount) +{ + + + Vector3f v0,v1,v2, tanu, tanv; + + Vector3f p0,p1,p2; + + Vector3f d1,d2; + + float uv[3][2],det,u,v,l1,l2; + + int i,j,k; + vector<Vector3f> sVector, tVector; + vector<Vector3f> avgS, avgT; + + avgS.resize (vertexCount); + avgT.resize (vertexCount); + + sVector.reserve (faceCount * 3); + tVector.reserve (faceCount * 3); + + for( i=0;i<vertexCount;i++ ) + + { + + sVector[i]=Vector3f(0.0, 0.0, 0.0); + + tVector[i]=Vector3f(0.0, 0.0, 0.0); + + } + + k=0; + for( i=0;i<faceCount;i++,k+=3 ) + + { + + v0=*((Vector3f *)&pPositions[indices[k]]); + + v1=*((Vector3f *)&pPositions[indices[k+1]])-v0; + + v2=*((Vector3f *)&pPositions[indices[k+2]])-v0; + + uv[0][0]=tex[indices[k]].x; + + uv[0][1]=tex[indices[k]].y; + + uv[1][0]=tex[indices[k+1]].x-uv[0][0]; + + uv[1][1]=tex[indices[k+1]].y-uv[0][1]; + + uv[2][0]=tex[indices[k+2]].x-uv[0][0]; + + uv[2][1]=tex[indices[k+2]].y-uv[0][1]; + + det=(uv[1][0]*uv[2][1])-(uv[2][0]*uv[1][1]); + + + + if (fabsf(det)<0.000001f){ + continue; + } + + u=0; v=0; + + u-=uv[0][0]; v-=uv[0][1]; + + p0=v0+v1*((u*uv[2][1]-uv[2][0]*v)/det)+v2*((uv[1][0]*v-u*uv[1][1])/det); + + + + u=1; v=0; + + u-=uv[0][0]; v-=uv[0][1]; + + p1=v0+v1*((u*uv[2][1]-uv[2][0]*v)/det)+v2*((uv[1][0]*v-u*uv[1][1])/det); + + + + u=0; v=1; + + u-=uv[0][0]; v-=uv[0][1]; + + p2=v0+v1*((u*uv[2][1]-uv[2][0]*v)/det)+v2*((uv[1][0]*v-u*uv[1][1])/det); + + + + d1=p2-p0; + + d2=p1-p0; + + l1=Magnitude(d1); + + l2=Magnitude(d2); + + d1*=1.0f/l1; + + d2*=1.0f/l2; + + + + j=indices[k]; + + sVector[j] += d1; + + tVector[j].x+=d2.x; tVector[j].y+=d2.y; tVector[j].z+=d2.z; + + + + j=indices[k+1]; + + sVector[j].x+=d1.x; sVector[j].y+=d1.y; sVector[j].z+=d1.z; + + tVector[j].x+=d2.x; tVector[j].y+=d2.y; tVector[j].z+=d2.z; + + + + j=indices[k+2]; + + sVector[j].x+=d1.x; sVector[j].y+=d1.y; sVector[j].z+=d1.z; + + tVector[j].x+=d2.x; tVector[j].y+=d2.y; tVector[j].z+=d2.z; + + } + + + + for( i=0;i<vertexCount;i++ ) + + { + + v0 = Vector3f(sVector[i].x,sVector[i].y,sVector[i].z); + v0 = NormalizeRobust(v0); + + v1 = Vector3f(tVector[i].x,tVector[i].y,tVector[i].z); + + Vector3f n(normals[i].x,normals[i].y,normals[i].z); + + if (SqrMagnitude(v1)<0.0001f) + { + v1 = Cross(v0,n); + } + v1 = NormalizeRobust(v1); + + + sVector[i]=v0; + + tVector[i].x=v1.x; + + tVector[i].y=v1.y; + + tVector[i].z=v1.z; + + avgOrthonormalBases[i].SetOrthoNormalBasis (tVector[i], sVector[i], n); + } + +} + diff --git a/Runtime/Geometry/TangentSpaceCalculation.h b/Runtime/Geometry/TangentSpaceCalculation.h new file mode 100644 index 0000000..2e925b0 --- /dev/null +++ b/Runtime/Geometry/TangentSpaceCalculation.h @@ -0,0 +1,13 @@ +#ifndef TANGENTSPACECALCULATION_H +#define TANGENTSPACECALCULATION_H + +class Vector3f; +class Vector2f; +class Matrix3x3f; + +void CreateTangentSpaceTangents (const Vector3f* pPositions, const Vector2f* tex, const Vector3f* normals, + Matrix3x3f* avgOrthonormalBases, int vertexCount, + const int* indices, int faceCount); + + +#endif diff --git a/Runtime/Geometry/TextureAtlas.cpp b/Runtime/Geometry/TextureAtlas.cpp new file mode 100644 index 0000000..6183bc3 --- /dev/null +++ b/Runtime/Geometry/TextureAtlas.cpp @@ -0,0 +1,588 @@ +#include "UnityPrefix.h" +#include "TextureAtlas.h" +#include "Runtime/Math/Rect.h" +#include "Runtime/Graphics/Texture2D.h" +#include "Runtime/Graphics/Image.h" +#include "Runtime/Math/Color.h" +#include "Runtime/Math/Vector2.h" +#include "Runtime/Utilities/BitUtility.h" +#include "Runtime/Shaders/GraphicsCaps.h" + +using namespace std; + +// Just implemented from this article: +// http://www.blackpawn.com/texts/lightmaps/default.html + +bool PackTextureAtlas( Texture2D* atlas, int atlasMaximumSize, int textureCount, Texture2D** textures, Rectf* outRects, int padding, int textureMode ); + +struct Node +{ + Node() : taken(false) { child[0] = NULL; child[1] = NULL; } + ~Node() { delete child[0]; delete child[1]; } + + void Reset() + { + delete child[0]; delete child[1]; + child[0] = NULL; child[1] = NULL; + taken = false; + } + + Node* Insert( float width, float height, float padding, bool use4PixelBoundaries ); + + Node* child[2]; + Rectf rect; + bool taken; +}; + + +Node* Node::Insert( float width, float height, float padding, bool use4PixelBoundaries ) +{ + // if we're not leaf, try inserting into children + if( child[0] ) + { + Node* newNode = child[0]->Insert( width, height, padding, use4PixelBoundaries ); + if( newNode ) + return newNode; + return child[1]->Insert( width, height, padding, use4PixelBoundaries ); + } + + // we are leaf + + if( taken ) + return NULL; // already taken + + // will it fit? + // 0.5 float error margin and we don't care about sub-texel overlaps anyway + if( width > rect.Width() - padding + 0.5f || height > rect.Height() - padding + 0.5f ) + return NULL; // won't fit + + // if this a perfect or nearly perfect fit, take it + float dw = rect.Width() - width; + float dh = rect.Height() - height; + if( dw <= padding*2 && dh <= padding*2 ) + { + taken = true; + return this; + } + if( use4PixelBoundaries && dw < 4 && dh < 4 ) + { + taken = true; + return this; + } + + // split the node + child[0] = new Node(); + child[1] = new Node(); + + // decide which way to split + if( dw > dh ) + { + // horizontal children + int split = int(width + padding); + if( use4PixelBoundaries ) + split = (split + 3) & (~3); + child[0]->rect = MinMaxRect( rect.x, rect.y, rect.x+width+padding, rect.GetBottom() ); + child[1]->rect = MinMaxRect( rect.x+split, rect.y, rect.GetRight(), rect.GetBottom() ); + } + else + { + // vertical children + int split = int(height + padding); + if( use4PixelBoundaries ) + split = (split + 3) & (~3); + child[0]->rect = MinMaxRect ( rect.x, rect.y, rect.GetRight(), rect.y+height+padding ); + child[1]->rect = MinMaxRect ( rect.x, rect.y+split, rect.GetRight(), rect.GetBottom() ); + } + + // insert into first child + return child[0]->Insert( width, height, padding, use4PixelBoundaries ); +} + +typedef std::pair<int,int> IntPair; +typedef std::vector<IntPair> TextureSizes; + +struct IndexSorter { + bool operator()( int a, int b ) const + { + return sizes[a].first * sizes[a].second > sizes[b].first * sizes[b].second; + } + IndexSorter( const TextureSizes& s ) : sizes(s) { } + const TextureSizes& sizes; +}; + +void PackAtlases (dynamic_array<Vector2f>& sizes, const int maxAtlasSize, const float padding, dynamic_array<Vector2f>& outOffsets, dynamic_array<int>& outIndices, int& atlasCount) +{ + const int count = sizes.size (); + const bool use4PixelBoundaries = false; + + dynamic_array<Node> atlases; + outOffsets.resize_uninitialized (count); + outIndices.resize_uninitialized (count); + + for (int i = 0; i < count; ++i) + { + Node* node = NULL; + int atlasIndex = -1; + while (!node) + { + atlasIndex++; + Vector2f& size = sizes[i]; + if (atlasIndex == atlases.size ()) + { + Node atlas; + atlas.rect.Set (0, 0, maxAtlasSize, maxAtlasSize); + atlases.push_back (atlas); + node = atlases[atlasIndex].Insert (size.x, size.y, padding, use4PixelBoundaries); + if (!node) + { + // We just tried inserting into an empty atlas. If that didn't succeed, we need to make the current rect smaller to fit maxAtlasSize + if (size.x > size.y) + { + size.y *= ((float)maxAtlasSize) / size.x; + size.x = maxAtlasSize; + } + else + { + size.x *= ((float)maxAtlasSize) / size.y; + size.y = maxAtlasSize; + } + node = atlases[atlasIndex].Insert (size.x, size.y, 0.0f, use4PixelBoundaries); + DebugAssert (node); + } + } + else + { + node = atlases[atlasIndex].Insert (size.x, size.y, padding, use4PixelBoundaries); + } + } + outOffsets[i].Set (node->rect.x, node->rect.y); + outIndices[i] = atlasIndex; + } + + atlasCount = atlases.size (); + + // deallocate all the trees + for (int i = 0; i < atlases.size (); ++i) + atlases[i].Reset (); +} + +bool PackTextureAtlasSimple( Texture2D* atlas, int atlasMaximumSize, int textureCount, Texture2D** textures, Rectf* outRects, int padding, bool upload, bool markNoLongerReadable ) +{ + atlasMaximumSize = min(gGraphicsCaps.maxTextureSize, atlasMaximumSize); + + // Cleanup the texture set. + // * Remove duplicate textures + // * Remove null textures + vector<int> remap; + remap.resize(textureCount); + vector<Texture2D*> uniqueTextures; + for (int i=0;i<textureCount;i++) + { + // Completely ignore null textures + if (textures[i] == NULL) + { + *outRects = Rectf (0,0,0,0); + remap[i] = -1; + continue; + } + + // Find duplicate texture and update remap + vector<Texture2D*> ::iterator found = find(uniqueTextures.begin(), uniqueTextures.end(), textures[i]); + if (found != uniqueTextures.end()) + { + remap[i] = distance(uniqueTextures.begin(), found); + } + else + { + remap[i] = uniqueTextures.size(); + uniqueTextures.push_back(textures[i]); + } + } + + if (!uniqueTextures.empty()) + { + vector<Rectf> uniqueRects; + uniqueRects.resize(uniqueTextures.size()); + + // Do the heavy lifting + if (!PackTextureAtlas(atlas, atlasMaximumSize, uniqueTextures.size(), &uniqueTextures[0], &uniqueRects[0], padding, upload ? 0 : Texture2D::kThreadedInitialize )) + return false; + + // Copy out rects from unique + for (int i=0;i<textureCount;i++) + { + if (remap[i] != -1) + outRects[i] = uniqueRects[remap[i]]; + } + } + + if (upload) + { + if (!IsAnyCompressedTextureFormat(atlas->GetTextureFormat())) + atlas->RebuildMipMap (); + + if( markNoLongerReadable ) + { + atlas->SetIsReadable(false); + atlas->SetIsUnreloadable(false); + } + + atlas->AwakeFromLoad(kDefaultAwakeFromLoad); + } + return true; +} + +enum PackingFormat { + kPackingUncompressed, + kPackingDXT1, + kPackingDXT5, +}; + + +bool PackTextureAtlas( Texture2D* atlas, int atlasMaximumSize, int textureCount, Texture2D** textures, Rectf* outRects, int padding, int textureOptions ) +{ + DebugAssertIf( !atlas || !textures || !outRects || textureCount <= 0 ); + const int kMinTextureSize = 4; + const int kMinAtlasSize = 8; + int i; + atlasMaximumSize = max (atlasMaximumSize, kMinAtlasSize); + + PackingFormat packFormat = kPackingDXT1; + bool packWithMipmaps = false; + bool someInputHasNoMipmaps = false; + + // Immediately decrease input texture sizes that are too large to fit; figure out + // result packing format and whether we'll have mipmaps. + TextureSizes textureSizes; + textureSizes.resize( textureCount ); + for( i = 0; i < textureCount; ++i ) + { + IntPair& size = textureSizes[i]; + size.first = textures[i]->GetDataWidth(); + size.second = textures[i]->GetDataHeight(); + while( size.first > atlasMaximumSize && size.first > kMinTextureSize ) + { + size.first /= 2; + packFormat = kPackingUncompressed; // we'll have to scale down, switch to uncompressed + } + while( size.second > atlasMaximumSize && size.second > kMinTextureSize ) + { + size.second /= 2; + packFormat = kPackingUncompressed; // we'll have to scale down, switch to uncompressed + } + + // Atlas format rules: + // Defaults to DXT1 + // If there is a DXT5 texture, pack to DXT5 (expand DXT1 alpha to opaque) + // If there is an uncompressed or DXT3 texture, pack to 32 bit uncompressed. + TextureFormat texFormat = textures[i]->GetTextureFormat(); + if (texFormat == kTexFormatDXT1 || texFormat == kTexFormatDXT5) + { + // Incoming texture is DXT1 or DXT5 + // If currently we are packing to DXT1 and incoming is DXT5, switch to that. + if( packFormat == kPackingDXT1 && texFormat == kTexFormatDXT5 ) + packFormat = kPackingDXT5; + } + else + { + // Incoming texture is anything else: pack to uncompressed + packFormat = kPackingUncompressed; + } + + // If any texture has mipmaps, then atlas will have them + if( textures[i]->HasMipMap() ) + packWithMipmaps = true; + else + someInputHasNoMipmaps = true; + } + + // If some input textures have mipmaps and some don't, then + // pack to uncompressed atlas. + if( packWithMipmaps && someInputHasNoMipmaps ) + packFormat = kPackingUncompressed; + + // Sort incoming textures by size; largest area first + std::vector<int> sortedIndices; + sortedIndices.resize( textureCount ); + for( i = 0; i < textureCount; ++i ) + { + sortedIndices[i] = i; + } + std::sort( sortedIndices.begin(), sortedIndices.end(), IndexSorter(textureSizes) ); + + // Calculate an initial lower bound estimate for the atlas width & height + int totalPixels = 0; + for( i = 0; i < textureCount; ++i ) + { + IntPair& size = textureSizes[i]; + totalPixels += size.first * size.second; + } + int atlasWidth = min<int>(NextPowerOfTwo(UInt32(Sqrt (totalPixels))), atlasMaximumSize); + int atlasHeight = min<int>(NextPowerOfTwo(totalPixels / atlasWidth), atlasMaximumSize); + // Do the packing of rectangles + bool packOk = true; + const int kMaxPackIterations = 100; + int packIterations = 0; + std::vector<Node*> textureNodes; + textureNodes.resize( textureCount ); + + // Create a tree to track occupied areas in the atlas + Node tree; + + do { + packOk = true; + tree.Reset(); + tree.rect = MinMaxRect<float> ( 0, 0, atlasWidth, atlasHeight ); + + bool use4PixelBoundaries = (packFormat != kPackingUncompressed); + + for( i = 0; i < textureCount; ++i ) + { + int texIndex = sortedIndices[i]; + DebugAssertIf( texIndex < 0 || texIndex >= textureCount ); + int texWidth = textureSizes[texIndex].first; + int texHeight = textureSizes[texIndex].second; + Node* node = tree.Insert( texWidth, texHeight, padding, use4PixelBoundaries ); + textureNodes[texIndex] = node; + if( !node ) + { + // texture does not fit; break out, reduce sizes and repack again + packOk = false; + break; + } + } + + // packing failed - decrease image sizes and try again + if( !packOk ) + { + // First we just increase the atlas size and see if we can fit all textures in. + if (atlasWidth != atlasMaximumSize || atlasHeight != atlasMaximumSize) + { + // Never increase beyond max size + if (atlasWidth == atlasMaximumSize) + atlasHeight *= 2; + else if (atlasHeight == atlasMaximumSize) + atlasWidth *= 2; + // increase the smaller of width/height + else if (atlasWidth < atlasHeight) + atlasWidth *= 2; + else + atlasHeight *= 2; + } + else + { + // TODO: the decreasing logic can be arbitrarily more complex. E.g. decrease the largest + // images first only, etc. + for( i = 0; i < textureCount; ++i ) + { + IntPair& size = textureSizes[i]; + if( size.first > kMinTextureSize && size.second > kMinTextureSize ) { + size.first = size.first * 3 / 4; + size.second = size.second * 3 / 4; + } + } + + // we'll scale images down, no DXT for ya + packFormat = kPackingUncompressed; + + // Only update pack iterations, for decreasing texture size + ++packIterations; + } + + AssertIf (atlasWidth > atlasMaximumSize); + AssertIf (atlasHeight > atlasMaximumSize); + } + } while( !packOk && packIterations < kMaxPackIterations ); + + if( !packOk ) + return false; + + + // Fill out UV rectangles for the input textures + for( i = 0; i < textureCount; ++i ) + { + int texIndex = sortedIndices[i]; + DebugAssertIf( texIndex < 0 || texIndex >= textureCount ); + int texWidth = textureSizes[texIndex].first; + int texHeight = textureSizes[texIndex].second; + const Node* node = textureNodes[texIndex]; + AssertIf( !node ); + + // Set the rectangle + outRects[texIndex] = MinMaxRect ( + node->rect.x/atlasWidth, + node->rect.y/atlasHeight, + (node->rect.x+texWidth)/atlasWidth, + (node->rect.y+texHeight)/atlasHeight ); + } + + + // Initialize atlas texture + TextureFormat atlasFormat; + if( packFormat == kPackingDXT1 ) + atlasFormat = kTexFormatDXT1; + else if( packFormat == kPackingDXT5 ) + atlasFormat = kTexFormatDXT5; + else + atlasFormat = kTexFormatARGB32; + + textureOptions |= packWithMipmaps ? (Texture2D::kMipmapMask) : (Texture2D::kNoMipmap); + atlas->InitTexture( atlasWidth, atlasHeight, atlasFormat, textureOptions ); + + // Packing into uncompressed texture atlas + if( packFormat == kPackingUncompressed ) + { + UInt8* atlasData = atlas->GetRawImageData(); + memset( atlasData, 0, atlas->GetRawImageData(1) - atlasData ); + + // Blit textures into final destinations + Image* decompressedImage = 0; + const int numAtlasMips = atlas->CountDataMipmaps(); + + for( i = 0; i < textureCount; ++i ) + { + int texIndex = sortedIndices[i]; + DebugAssertIf( texIndex < 0 || texIndex >= textureCount ); + + int texWidth = textureSizes[texIndex].first; + int texHeight = textureSizes[texIndex].second; + const Node* node = textureNodes[texIndex]; + AssertIf( !node ); + + int destCoordX = node->rect.x; + int destCoordY = node->rect.y; + int destWidth = std::min(texWidth, std::max(1, (int)node->rect.width - padding)); + int destHeight = std::min(texHeight, std::max(1, (int)node->rect.height - padding)); + + int atlasMipWidth = atlasWidth; + int atlasMipHeight = atlasHeight; + // copy all mips to atlas + for ( int mip=0, numMips=std::min( textures[texIndex]->CountDataMipmaps(), numAtlasMips ); mip!=numMips; ++mip ) + { + ImageReference atlasImgRef; + atlas->GetWriteImageReference ( &atlasImgRef, 0, mip ); + // get texture rect in atlas for current source texture + ImageReference destRect = atlasImgRef.ClipImage( destCoordX, destCoordY, destWidth, destHeight ); + + ImageReference srcMip; + if ( textures[texIndex]->GetWriteImageReference( &srcMip, 0, mip ) ) + { + ImageReference::BlitMode blit_mode = ImageReference::BLIT_BILINEAR_SCALE; + if ( destWidth==srcMip.GetWidth() && destHeight==srcMip.GetHeight() ) + blit_mode = ImageReference::BLIT_COPY; + destRect.BlitImage( srcMip, blit_mode ); + } + else + { + if( !decompressedImage ) + decompressedImage = new Image( texWidth, texHeight, kTexFormatRGBA32 ); + else + decompressedImage->SetImage( texWidth, texHeight, kTexFormatRGBA32, false ); + + textures[texIndex]->ExtractImage( decompressedImage, 0 ); + ImageReference::BlitMode blit_mode = ImageReference::BLIT_BILINEAR_SCALE; + if ( destWidth==decompressedImage->GetWidth() && destHeight==decompressedImage->GetHeight() ) + blit_mode = ImageReference::BLIT_COPY; + destRect.BlitImage( *decompressedImage, blit_mode ); + } + + // Go to next mip level + destCoordX /= 2; + destCoordY /= 2; + destWidth /= 2; + destHeight /= 2; + atlasMipWidth = std::max( atlasMipWidth/2, 1 ); + atlasMipHeight = std::max( atlasMipHeight/2, 1 ); + texWidth = std::max( texWidth/2, 1 ); + texHeight = std::max( texHeight/2, 1 ); + } + } + delete decompressedImage; + } + // Packing into compressed texture atlas + else + { + UInt8* atlasData = atlas->GetRawImageData(); + bool atlasDXT1 = (atlasFormat==kTexFormatDXT1); + int blockBytes = atlasDXT1 ? 8 : 16; + memset( atlasData, 0, atlas->GetRawImageData(1) - atlasData ); + + // Blit textures into final destinations + for( i = 0; i < textureCount; ++i ) + { + int texIndex = sortedIndices[i]; + DebugAssertIf( texIndex < 0 || texIndex >= textureCount ); + int texWidth = textureSizes[texIndex].first; + int texHeight = textureSizes[texIndex].second; + Texture2D* src = textures[texIndex]; + int mipCount = std::min( src->CountDataMipmaps(), atlas->CountDataMipmaps() ); + AssertIf( texWidth != src->GetDataWidth() || texHeight != src->GetDataHeight() ); + const Node* node = textureNodes[texIndex]; + AssertIf( !node ); + int destCoordX = int(node->rect.x), destCoordY = int(node->rect.y); + int destWidth = int(node->rect.width), destHeight = int(node->rect.height); + AssertIf( (destCoordX & 3) != 0 || (destCoordY & 3) != 0 ); + + UInt8* atlasMipData = atlasData; + int atlasMipWidth = atlasWidth; + int atlasMipHeight = atlasHeight; + const UInt8* srcPointer = src->GetRawImageData(); + + if(srcPointer) + { + + // blit mip levels while we have them + for( int mip = 0; mip < mipCount; ++mip ) + { + // Get pointer where should we blit texture into + int destBlockX = destCoordX / 4; + int destBlockY = destCoordY / 4; + UInt8* destPointer = atlasMipData + (destBlockY * atlasMipWidth/4 + destBlockX) * blockBytes; + + TextureFormat srcFormat = src->GetTextureFormat(); + AssertIf( !IsCompressedDXTTextureFormat(srcFormat) ); + if( srcFormat == atlasFormat ) + { + BlitCopyCompressedImage( srcFormat, srcPointer, + texWidth, texHeight, + destPointer, atlasMipWidth, atlasMipHeight, false ); + } + else if( srcFormat == kTexFormatDXT1 && atlasFormat == kTexFormatDXT5 ) + { + BlitCopyCompressedDXT1ToDXT5( srcPointer, + texWidth, texHeight, + destPointer, atlasMipWidth, atlasMipHeight ); + } + else + { + AssertString( "Unsupported format in compressed texture atlas" ); + } + + // Go to next mip level + srcPointer += CalculateImageSize( texWidth, texHeight, srcFormat ); + atlasMipData += CalculateImageSize( atlasMipWidth, atlasMipHeight, atlasFormat ); + destCoordX /= 2; + destCoordY /= 2; + destWidth /= 2; + destHeight /= 2; + atlasMipWidth = std::max( atlasMipWidth / 2, 4 ); + atlasMipHeight = std::max( atlasMipHeight / 2, 4 ); + texWidth = std::max( texWidth / 2, 4 ); + texHeight = std::max( texHeight / 2, 4 ); + + // Stop if we begin to straddle DXT block boundaries. + if( (destCoordX & 3) != 0 || (destCoordY & 3) != 0 ) + break; + + // Stop if we don't fit into our initial area anymore. + if( destWidth < 4 || destHeight < 4 ) + break; + } + } + else + ErrorStringMsg ("Could not read texture data for texture '%s'. Make sure that Read/Write access is enabled in the texture importer advanced settings\n", src->GetName()); + } + } + + return true; +} diff --git a/Runtime/Geometry/TextureAtlas.h b/Runtime/Geometry/TextureAtlas.h new file mode 100644 index 0000000..132f8ef --- /dev/null +++ b/Runtime/Geometry/TextureAtlas.h @@ -0,0 +1,12 @@ +#pragma once + +#include "Runtime/Utilities/dynamic_array.h" +#include "Runtime/Math/Vector2.h" +#include "Runtime/Math/Rect.h" +#include "Runtime/Modules/ExportModules.h" + +class Texture2D; + +bool EXPORT_COREMODULE PackTextureAtlasSimple( Texture2D* atlas, int atlasMaximumSize, int textureCount, Texture2D** textures, Rectf* outRects, int padding, bool upload, bool markNoLongerReadable ); +void PackAtlases (dynamic_array<Vector2f>& sizes, const int maxAtlasSize, const float padding, dynamic_array<Vector2f>& outOffsets, dynamic_array<int>& outIndices, int& atlasCount); + diff --git a/Runtime/Geometry/TriTriIntersect.cpp b/Runtime/Geometry/TriTriIntersect.cpp new file mode 100644 index 0000000..03dfdc4 --- /dev/null +++ b/Runtime/Geometry/TriTriIntersect.cpp @@ -0,0 +1,719 @@ +#include "UnityPrefix.h" +/* Triangle/triangle intersection test routine, + * by Tomas Moller, 1997. + * See article "A Fast Triangle-Triangle Intersection Test", + * Journal of Graphics Tools, 2(2), 1997 + * updated: 2001-06-20 (added line of intersection) + * + * int tri_tri_intersect(float V0[3],float V1[3],float V2[3], + * float U0[3],float U1[3],float U2[3]) + * + * parameters: vertices of triangle 1: V0,V1,V2 + * vertices of triangle 2: U0,U1,U2 + * result : returns 1 if the triangles intersect, otherwise 0 + * + * Here is a version withouts divisions (a little faster) + * int NoDivTriTriIsect(float V0[3],float V1[3],float V2[3], + * float U0[3],float U1[3],float U2[3]); + * + * This version computes the line of intersection as well (if they are not coplanar): + * int tri_tri_intersect_with_isectline(float V0[3],float V1[3],float V2[3], + * float U0[3],float U1[3],float U2[3],int *coplanar, + * float isectpt1[3],float isectpt2[3]); + * coplanar returns whether the tris are coplanar + * isectpt1, isectpt2 are the endpoints of the line of intersection + */ + +#include <math.h> +static int coplanar_tri_tri(float N[3],float V0[3],float V1[3],float V2[3], + float U0[3],float U1[3],float U2[3]); +static int tri_tri_intersect(float V0[3],float V1[3],float V2[3], + float U0[3],float U1[3],float U2[3]); +static int NoDivTriTriIsect(float V0[3],float V1[3],float V2[3], + float U0[3],float U1[3],float U2[3]); +int tri_tri_intersect_with_isectline(float V0[3],float V1[3],float V2[3], + float U0[3],float U1[3],float U2[3],int *coplanar, + float isectpt1[3],float isectpt2[3]); + +#define FABS(x) ((float)fabs(x)) /* implement as is fastest on your machine */ + +/* if USE_EPSILON_TEST is true then we do a check: + if |dv|<EPSILON then dv=0.0; + else no check is done (which is less robust) +*/ +#ifndef TRUE +#define TRUE 1 +#endif + +#define USE_EPSILON_TEST TRUE +#define EPSILON 0.000001 + + +/* some macros */ +#define CROSS(dest,v1,v2) \ + dest[0]=v1[1]*v2[2]-v1[2]*v2[1]; \ + dest[1]=v1[2]*v2[0]-v1[0]*v2[2]; \ + dest[2]=v1[0]*v2[1]-v1[1]*v2[0]; + +#define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]) + +#define SUB(dest,v1,v2) dest[0]=v1[0]-v2[0]; dest[1]=v1[1]-v2[1]; dest[2]=v1[2]-v2[2]; + +#define ADD(dest,v1,v2) dest[0]=v1[0]+v2[0]; dest[1]=v1[1]+v2[1]; dest[2]=v1[2]+v2[2]; + +#define MULT(dest,v,factor) dest[0]=factor*v[0]; dest[1]=factor*v[1]; dest[2]=factor*v[2]; + +#define SET(dest,src) dest[0]=src[0]; dest[1]=src[1]; dest[2]=src[2]; + +/* sort so that a<=b */ +#define SORT(a,b) \ + if(a>b) \ + { \ + float c; \ + c=a; \ + a=b; \ + b=c; \ + } + +#define ISECT(VV0,VV1,VV2,D0,D1,D2,isect0,isect1) \ + isect0=VV0+(VV1-VV0)*D0/(D0-D1); \ + isect1=VV0+(VV2-VV0)*D0/(D0-D2); + + +#define COMPUTE_INTERVALS(VV0,VV1,VV2,D0,D1,D2,D0D1,D0D2,isect0,isect1) \ + if(D0D1>0.0f) \ + { \ + /* here we know that D0D2<=0.0 */ \ + /* that is D0, D1 are on the same side, D2 on the other or on the plane */ \ + ISECT(VV2,VV0,VV1,D2,D0,D1,isect0,isect1); \ + } \ + else if(D0D2>0.0f) \ + { \ + /* here we know that d0d1<=0.0 */ \ + ISECT(VV1,VV0,VV2,D1,D0,D2,isect0,isect1); \ + } \ + else if(D1*D2>0.0f || D0!=0.0f) \ + { \ + /* here we know that d0d1<=0.0 or that D0!=0.0 */ \ + ISECT(VV0,VV1,VV2,D0,D1,D2,isect0,isect1); \ + } \ + else if(D1!=0.0f) \ + { \ + ISECT(VV1,VV0,VV2,D1,D0,D2,isect0,isect1); \ + } \ + else if(D2!=0.0f) \ + { \ + ISECT(VV2,VV0,VV1,D2,D0,D1,isect0,isect1); \ + } \ + else \ + { \ + /* triangles are coplanar */ \ + return coplanar_tri_tri(N1,V0,V1,V2,U0,U1,U2); \ + } + + + +/* this edge to edge test is based on Franlin Antonio's gem: + "Faster Line Segment Intersection", in Graphics Gems III, + pp. 199-202 */ +#define EDGE_EDGE_TEST(V0,U0,U1) \ + Bx=U0[i0]-U1[i0]; \ + By=U0[i1]-U1[i1]; \ + Cx=V0[i0]-U0[i0]; \ + Cy=V0[i1]-U0[i1]; \ + f=Ay*Bx-Ax*By; \ + d=By*Cx-Bx*Cy; \ + if((f>0 && d>=0 && d<=f) || (f<0 && d<=0 && d>=f)) \ + { \ + e=Ax*Cy-Ay*Cx; \ + if(f>0) \ + { \ + if(e>=0 && e<=f) return 1; \ + } \ + else \ + { \ + if(e<=0 && e>=f) return 1; \ + } \ + } + +#define EDGE_AGAINST_TRI_EDGES(V0,V1,U0,U1,U2) \ +{ \ + float Ax,Ay,Bx,By,Cx,Cy,e,d,f; \ + Ax=V1[i0]-V0[i0]; \ + Ay=V1[i1]-V0[i1]; \ + /* test edge U0,U1 against V0,V1 */ \ + EDGE_EDGE_TEST(V0,U0,U1); \ + /* test edge U1,U2 against V0,V1 */ \ + EDGE_EDGE_TEST(V0,U1,U2); \ + /* test edge U2,U1 against V0,V1 */ \ + EDGE_EDGE_TEST(V0,U2,U0); \ +} + +#define POINT_IN_TRI(V0,U0,U1,U2) \ +{ \ + float a,b,c,d0,d1,d2; \ + /* is T1 completly inside T2? */ \ + /* check if V0 is inside tri(U0,U1,U2) */ \ + a=U1[i1]-U0[i1]; \ + b=-(U1[i0]-U0[i0]); \ + c=-a*U0[i0]-b*U0[i1]; \ + d0=a*V0[i0]+b*V0[i1]+c; \ + \ + a=U2[i1]-U1[i1]; \ + b=-(U2[i0]-U1[i0]); \ + c=-a*U1[i0]-b*U1[i1]; \ + d1=a*V0[i0]+b*V0[i1]+c; \ + \ + a=U0[i1]-U2[i1]; \ + b=-(U0[i0]-U2[i0]); \ + c=-a*U2[i0]-b*U2[i1]; \ + d2=a*V0[i0]+b*V0[i1]+c; \ + if(d0*d1>0.0) \ + { \ + if(d0*d2>0.0) return 1; \ + } \ +} + +static int coplanar_tri_tri(float N[3],float V0[3],float V1[3],float V2[3], + float U0[3],float U1[3],float U2[3]) +{ + float A[3]; + short i0,i1; + /* first project onto an axis-aligned plane, that maximizes the area */ + /* of the triangles, compute indices: i0,i1. */ + A[0]=fabs(N[0]); + A[1]=fabs(N[1]); + A[2]=fabs(N[2]); + if(A[0]>A[1]) + { + if(A[0]>A[2]) + { + i0=1; /* A[0] is greatest */ + i1=2; + } + else + { + i0=0; /* A[2] is greatest */ + i1=1; + } + } + else /* A[0]<=A[1] */ + { + if(A[2]>A[1]) + { + i0=0; /* A[2] is greatest */ + i1=1; + } + else + { + i0=0; /* A[1] is greatest */ + i1=2; + } + } + + /* test all edges of triangle 1 against the edges of triangle 2 */ + EDGE_AGAINST_TRI_EDGES(V0,V1,U0,U1,U2); + EDGE_AGAINST_TRI_EDGES(V1,V2,U0,U1,U2); + EDGE_AGAINST_TRI_EDGES(V2,V0,U0,U1,U2); + + /* finally, test if tri1 is totally contained in tri2 or vice versa */ + POINT_IN_TRI(V0,U0,U1,U2); + POINT_IN_TRI(U0,V0,V1,V2); + + return 0; +} + + +static int tri_tri_intersect(float V0[3],float V1[3],float V2[3], + float U0[3],float U1[3],float U2[3]) +{ + float E1[3],E2[3]; + float N1[3],N2[3],d1,d2; + float du0,du1,du2,dv0,dv1,dv2; + float D[3]; + float isect1[2], isect2[2]; + float du0du1,du0du2,dv0dv1,dv0dv2; + short index; + float vp0,vp1,vp2; + float up0,up1,up2; + float b,c,max; + + /* compute plane equation of triangle(V0,V1,V2) */ + SUB(E1,V1,V0); + SUB(E2,V2,V0); + CROSS(N1,E1,E2); + d1=-DOT(N1,V0); + /* plane equation 1: N1.X+d1=0 */ + + /* put U0,U1,U2 into plane equation 1 to compute signed distances to the plane*/ + du0=DOT(N1,U0)+d1; + du1=DOT(N1,U1)+d1; + du2=DOT(N1,U2)+d1; + + /* coplanarity robustness check */ +#if USE_EPSILON_TEST==TRUE + if(fabs(du0)<EPSILON) du0=0.0; + if(fabs(du1)<EPSILON) du1=0.0; + if(fabs(du2)<EPSILON) du2=0.0; +#endif + du0du1=du0*du1; + du0du2=du0*du2; + + if(du0du1>0.0f && du0du2>0.0f) /* same sign on all of them + not equal 0 ? */ + return 0; /* no intersection occurs */ + + /* compute plane of triangle (U0,U1,U2) */ + SUB(E1,U1,U0); + SUB(E2,U2,U0); + CROSS(N2,E1,E2); + d2=-DOT(N2,U0); + /* plane equation 2: N2.X+d2=0 */ + + /* put V0,V1,V2 into plane equation 2 */ + dv0=DOT(N2,V0)+d2; + dv1=DOT(N2,V1)+d2; + dv2=DOT(N2,V2)+d2; + +#if USE_EPSILON_TEST==TRUE + if(fabs(dv0)<EPSILON) dv0=0.0; + if(fabs(dv1)<EPSILON) dv1=0.0; + if(fabs(dv2)<EPSILON) dv2=0.0; +#endif + + dv0dv1=dv0*dv1; + dv0dv2=dv0*dv2; + + if(dv0dv1>0.0f && dv0dv2>0.0f) /* same sign on all of them + not equal 0 ? */ + return 0; /* no intersection occurs */ + + /* compute direction of intersection line */ + CROSS(D,N1,N2); + + /* compute and index to the largest component of D */ + max=fabs(D[0]); + index=0; + b=fabs(D[1]); + c=fabs(D[2]); + if(b>max) max=b,index=1; + if(c>max) max=c,index=2; + + /* this is the simplified projection onto L*/ + vp0=V0[index]; + vp1=V1[index]; + vp2=V2[index]; + + up0=U0[index]; + up1=U1[index]; + up2=U2[index]; + + /* compute interval for triangle 1 */ + COMPUTE_INTERVALS(vp0,vp1,vp2,dv0,dv1,dv2,dv0dv1,dv0dv2,isect1[0],isect1[1]); + + /* compute interval for triangle 2 */ + COMPUTE_INTERVALS(up0,up1,up2,du0,du1,du2,du0du1,du0du2,isect2[0],isect2[1]); + + SORT(isect1[0],isect1[1]); + SORT(isect2[0],isect2[1]); + + if(isect1[1]<isect2[0] || isect2[1]<isect1[0]) return 0; + return 1; +} + + +#define NEWCOMPUTE_INTERVALS(VV0,VV1,VV2,D0,D1,D2,D0D1,D0D2,A,B,C,X0,X1) \ +{ \ + if(D0D1>0.0f) \ + { \ + /* here we know that D0D2<=0.0 */ \ + /* that is D0, D1 are on the same side, D2 on the other or on the plane */ \ + A=VV2; B=(VV0-VV2)*D2; C=(VV1-VV2)*D2; X0=D2-D0; X1=D2-D1; \ + } \ + else if(D0D2>0.0f)\ + { \ + /* here we know that d0d1<=0.0 */ \ + A=VV1; B=(VV0-VV1)*D1; C=(VV2-VV1)*D1; X0=D1-D0; X1=D1-D2; \ + } \ + else if(D1*D2>0.0f || D0!=0.0f) \ + { \ + /* here we know that d0d1<=0.0 or that D0!=0.0 */ \ + A=VV0; B=(VV1-VV0)*D0; C=(VV2-VV0)*D0; X0=D0-D1; X1=D0-D2; \ + } \ + else if(D1!=0.0f) \ + { \ + A=VV1; B=(VV0-VV1)*D1; C=(VV2-VV1)*D1; X0=D1-D0; X1=D1-D2; \ + } \ + else if(D2!=0.0f) \ + { \ + A=VV2; B=(VV0-VV2)*D2; C=(VV1-VV2)*D2; X0=D2-D0; X1=D2-D1; \ + } \ + else \ + { \ + /* triangles are coplanar */ \ + return coplanar_tri_tri(N1,V0,V1,V2,U0,U1,U2); \ + } \ +} + + + +static int NoDivTriTriIsect(float V0[3],float V1[3],float V2[3], + float U0[3],float U1[3],float U2[3]) +{ + float E1[3],E2[3]; + float N1[3],N2[3],d1,d2; + float du0,du1,du2,dv0,dv1,dv2; + float D[3]; + float isect1[2], isect2[2]; + float du0du1,du0du2,dv0dv1,dv0dv2; + short index; + float vp0,vp1,vp2; + float up0,up1,up2; + float bb,cc,max; + float a,b,c,x0,x1; + float d,e,f,y0,y1; + float xx,yy,xxyy,tmp; + + /* compute plane equation of triangle(V0,V1,V2) */ + SUB(E1,V1,V0); + SUB(E2,V2,V0); + CROSS(N1,E1,E2); + d1=-DOT(N1,V0); + /* plane equation 1: N1.X+d1=0 */ + + /* put U0,U1,U2 into plane equation 1 to compute signed distances to the plane*/ + du0=DOT(N1,U0)+d1; + du1=DOT(N1,U1)+d1; + du2=DOT(N1,U2)+d1; + + /* coplanarity robustness check */ +#if USE_EPSILON_TEST==TRUE + if(FABS(du0)<EPSILON) du0=0.0; + if(FABS(du1)<EPSILON) du1=0.0; + if(FABS(du2)<EPSILON) du2=0.0; +#endif + du0du1=du0*du1; + du0du2=du0*du2; + + if(du0du1>0.0f && du0du2>0.0f) /* same sign on all of them + not equal 0 ? */ + return 0; /* no intersection occurs */ + + /* compute plane of triangle (U0,U1,U2) */ + SUB(E1,U1,U0); + SUB(E2,U2,U0); + CROSS(N2,E1,E2); + d2=-DOT(N2,U0); + /* plane equation 2: N2.X+d2=0 */ + + /* put V0,V1,V2 into plane equation 2 */ + dv0=DOT(N2,V0)+d2; + dv1=DOT(N2,V1)+d2; + dv2=DOT(N2,V2)+d2; + +#if USE_EPSILON_TEST==TRUE + if(FABS(dv0)<EPSILON) dv0=0.0; + if(FABS(dv1)<EPSILON) dv1=0.0; + if(FABS(dv2)<EPSILON) dv2=0.0; +#endif + + dv0dv1=dv0*dv1; + dv0dv2=dv0*dv2; + + if(dv0dv1>0.0f && dv0dv2>0.0f) /* same sign on all of them + not equal 0 ? */ + return 0; /* no intersection occurs */ + + /* compute direction of intersection line */ + CROSS(D,N1,N2); + + /* compute and index to the largest component of D */ + max=(float)FABS(D[0]); + index=0; + bb=(float)FABS(D[1]); + cc=(float)FABS(D[2]); + if(bb>max) max=bb,index=1; + if(cc>max) max=cc,index=2; + + /* this is the simplified projection onto L*/ + vp0=V0[index]; + vp1=V1[index]; + vp2=V2[index]; + + up0=U0[index]; + up1=U1[index]; + up2=U2[index]; + + /* compute interval for triangle 1 */ + NEWCOMPUTE_INTERVALS(vp0,vp1,vp2,dv0,dv1,dv2,dv0dv1,dv0dv2,a,b,c,x0,x1); + + /* compute interval for triangle 2 */ + NEWCOMPUTE_INTERVALS(up0,up1,up2,du0,du1,du2,du0du1,du0du2,d,e,f,y0,y1); + + xx=x0*x1; + yy=y0*y1; + xxyy=xx*yy; + + tmp=a*xxyy; + isect1[0]=tmp+b*x1*yy; + isect1[1]=tmp+c*x0*yy; + + tmp=d*xxyy; + isect2[0]=tmp+e*xx*y1; + isect2[1]=tmp+f*xx*y0; + + SORT(isect1[0],isect1[1]); + SORT(isect2[0],isect2[1]); + + if(isect1[1]<isect2[0] || isect2[1]<isect1[0]) return 0; + return 1; +} + +/* sort so that a<=b */ +#define SORT2(a,b,smallest) \ + if(a>b) \ + { \ + float c; \ + c=a; \ + a=b; \ + b=c; \ + smallest=1; \ + } \ + else smallest=0; + + +inline void isect2(float VTX0[3],float VTX1[3],float VTX2[3],float VV0,float VV1,float VV2, + float D0,float D1,float D2,float *isect0,float *isect1,float isectpoint0[3],float isectpoint1[3]) +{ + float tmp=D0/(D0-D1); + float diff[3]; + *isect0=VV0+(VV1-VV0)*tmp; + SUB(diff,VTX1,VTX0); + MULT(diff,diff,tmp); + ADD(isectpoint0,diff,VTX0); + tmp=D0/(D0-D2); + *isect1=VV0+(VV2-VV0)*tmp; + SUB(diff,VTX2,VTX0); + MULT(diff,diff,tmp); + ADD(isectpoint1,VTX0,diff); +} + + +#if 0 +#define ISECT2(VTX0,VTX1,VTX2,VV0,VV1,VV2,D0,D1,D2,isect0,isect1,isectpoint0,isectpoint1) \ + tmp=D0/(D0-D1); \ + isect0=VV0+(VV1-VV0)*tmp; \ + SUB(diff,VTX1,VTX0); \ + MULT(diff,diff,tmp); \ + ADD(isectpoint0,diff,VTX0); \ + tmp=D0/(D0-D2); +/* isect1=VV0+(VV2-VV0)*tmp; \ */ +/* SUB(diff,VTX2,VTX0); \ */ +/* MULT(diff,diff,tmp); \ */ +/* ADD(isectpoint1,VTX0,diff); */ +#endif + +inline int compute_intervals_isectline(float VERT0[3],float VERT1[3],float VERT2[3], + float VV0,float VV1,float VV2,float D0,float D1,float D2, + float D0D1,float D0D2,float *isect0,float *isect1, + float isectpoint0[3],float isectpoint1[3]) +{ + if(D0D1>0.0f) + { + /* here we know that D0D2<=0.0 */ + /* that is D0, D1 are on the same side, D2 on the other or on the plane */ + isect2(VERT2,VERT0,VERT1,VV2,VV0,VV1,D2,D0,D1,isect0,isect1,isectpoint0,isectpoint1); + } + else if(D0D2>0.0f) + { + /* here we know that d0d1<=0.0 */ + isect2(VERT1,VERT0,VERT2,VV1,VV0,VV2,D1,D0,D2,isect0,isect1,isectpoint0,isectpoint1); + } + else if(D1*D2>0.0f || D0!=0.0f) + { + /* here we know that d0d1<=0.0 or that D0!=0.0 */ + isect2(VERT0,VERT1,VERT2,VV0,VV1,VV2,D0,D1,D2,isect0,isect1,isectpoint0,isectpoint1); + } + else if(D1!=0.0f) + { + isect2(VERT1,VERT0,VERT2,VV1,VV0,VV2,D1,D0,D2,isect0,isect1,isectpoint0,isectpoint1); + } + else if(D2!=0.0f) + { + isect2(VERT2,VERT0,VERT1,VV2,VV0,VV1,D2,D0,D1,isect0,isect1,isectpoint0,isectpoint1); + } + else + { + /* triangles are coplanar */ + return 1; + } + return 0; +} + +#define COMPUTE_INTERVALS_ISECTLINE(VERT0,VERT1,VERT2,VV0,VV1,VV2,D0,D1,D2,D0D1,D0D2,isect0,isect1,isectpoint0,isectpoint1) \ + if(D0D1>0.0f) \ + { \ + /* here we know that D0D2<=0.0 */ \ + /* that is D0, D1 are on the same side, D2 on the other or on the plane */ \ + isect2(VERT2,VERT0,VERT1,VV2,VV0,VV1,D2,D0,D1,&isect0,&isect1,isectpoint0,isectpoint1); \ + } +#if 0 + else if(D0D2>0.0f) \ + { \ + /* here we know that d0d1<=0.0 */ \ + isect2(VERT1,VERT0,VERT2,VV1,VV0,VV2,D1,D0,D2,&isect0,&isect1,isectpoint0,isectpoint1); \ + } \ + else if(D1*D2>0.0f || D0!=0.0f) \ + { \ + /* here we know that d0d1<=0.0 or that D0!=0.0 */ \ + isect2(VERT0,VERT1,VERT2,VV0,VV1,VV2,D0,D1,D2,&isect0,&isect1,isectpoint0,isectpoint1); \ + } \ + else if(D1!=0.0f) \ + { \ + isect2(VERT1,VERT0,VERT2,VV1,VV0,VV2,D1,D0,D2,&isect0,&isect1,isectpoint0,isectpoint1); \ + } \ + else if(D2!=0.0f) \ + { \ + isect2(VERT2,VERT0,VERT1,VV2,VV0,VV1,D2,D0,D1,&isect0,&isect1,isectpoint0,isectpoint1); \ + } \ + else \ + { \ + /* triangles are coplanar */ \ + coplanar=1; \ + return coplanar_tri_tri(N1,V0,V1,V2,U0,U1,U2); \ + } +#endif + +int tri_tri_intersect_with_isectline(float V0[3],float V1[3],float V2[3], + float U0[3],float U1[3],float U2[3],int *coplanar, + float isectpt1[3],float isectpt2[3]) +{ + float E1[3],E2[3]; + float N1[3],N2[3],d1,d2; + float du0,du1,du2,dv0,dv1,dv2; + float D[3]; + float isect1[2], isect2[2] = {.0f, .0f}; + float isectpointA1[3],isectpointA2[3]; + float isectpointB1[3] = {.0f, .0f, .0f},isectpointB2[3] = {.0f, .0f, .0f}; + float du0du1,du0du2,dv0dv1,dv0dv2; + short index; + float vp0,vp1,vp2; + float up0,up1,up2; + float b,c,max; + int smallest1,smallest2; + + /* compute plane equation of triangle(V0,V1,V2) */ + SUB(E1,V1,V0); + SUB(E2,V2,V0); + CROSS(N1,E1,E2); + d1=-DOT(N1,V0); + /* plane equation 1: N1.X+d1=0 */ + + /* put U0,U1,U2 into plane equation 1 to compute signed distances to the plane*/ + du0=DOT(N1,U0)+d1; + du1=DOT(N1,U1)+d1; + du2=DOT(N1,U2)+d1; + + /* coplanarity robustness check */ +#if USE_EPSILON_TEST==TRUE + if(fabs(du0)<EPSILON) du0=0.0; + if(fabs(du1)<EPSILON) du1=0.0; + if(fabs(du2)<EPSILON) du2=0.0; +#endif + du0du1=du0*du1; + du0du2=du0*du2; + + if(du0du1>0.0f && du0du2>0.0f) /* same sign on all of them + not equal 0 ? */ + return 0; /* no intersection occurs */ + + /* compute plane of triangle (U0,U1,U2) */ + SUB(E1,U1,U0); + SUB(E2,U2,U0); + CROSS(N2,E1,E2); + d2=-DOT(N2,U0); + /* plane equation 2: N2.X+d2=0 */ + + /* put V0,V1,V2 into plane equation 2 */ + dv0=DOT(N2,V0)+d2; + dv1=DOT(N2,V1)+d2; + dv2=DOT(N2,V2)+d2; + +#if USE_EPSILON_TEST==TRUE + if(fabs(dv0)<EPSILON) dv0=0.0; + if(fabs(dv1)<EPSILON) dv1=0.0; + if(fabs(dv2)<EPSILON) dv2=0.0; +#endif + + dv0dv1=dv0*dv1; + dv0dv2=dv0*dv2; + + if(dv0dv1>0.0f && dv0dv2>0.0f) /* same sign on all of them + not equal 0 ? */ + return 0; /* no intersection occurs */ + + /* compute direction of intersection line */ + CROSS(D,N1,N2); + + /* compute and index to the largest component of D */ + max=fabs(D[0]); + index=0; + b=fabs(D[1]); + c=fabs(D[2]); + if(b>max) max=b,index=1; + if(c>max) max=c,index=2; + + /* this is the simplified projection onto L*/ + vp0=V0[index]; + vp1=V1[index]; + vp2=V2[index]; + + up0=U0[index]; + up1=U1[index]; + up2=U2[index]; + + /* compute interval for triangle 1 */ + *coplanar=compute_intervals_isectline(V0,V1,V2,vp0,vp1,vp2,dv0,dv1,dv2, + dv0dv1,dv0dv2,&isect1[0],&isect1[1],isectpointA1,isectpointA2); + if(*coplanar) return coplanar_tri_tri(N1,V0,V1,V2,U0,U1,U2); + + + /* compute interval for triangle 2 */ + compute_intervals_isectline(U0,U1,U2,up0,up1,up2,du0,du1,du2, + du0du1,du0du2,&isect2[0],&isect2[1],isectpointB1,isectpointB2); + + SORT2(isect1[0],isect1[1],smallest1); + SORT2(isect2[0],isect2[1],smallest2); + + if(isect1[1]<isect2[0] || isect2[1]<isect1[0]) return 0; + + /* at this point, we know that the triangles intersect */ + + if(isect2[0]<isect1[0]) + { + if(smallest1==0) { SET(isectpt1,isectpointA1); } + else { SET(isectpt1,isectpointA2); } + + if(isect2[1]<isect1[1]) + { + if(smallest2==0) { SET(isectpt2,isectpointB2); } + else { SET(isectpt2,isectpointB1); } + } + else + { + if(smallest1==0) { SET(isectpt2,isectpointA2); } + else { SET(isectpt2,isectpointA1); } + } + } + else + { + if(smallest2==0) { SET(isectpt1,isectpointB1); } + else { SET(isectpt1,isectpointB2); } + + if(isect2[1]>isect1[1]) + { + if(smallest1==0) { SET(isectpt2,isectpointA2); } + else { SET(isectpt2,isectpointA1); } + } + else + { + if(smallest2==0) { SET(isectpt2,isectpointB2); } + else { SET(isectpt2,isectpointB1); } + } + } + return 1; +} diff --git a/Runtime/Geometry/TriTriIntersect.h b/Runtime/Geometry/TriTriIntersect.h new file mode 100644 index 0000000..ef3ecb4 --- /dev/null +++ b/Runtime/Geometry/TriTriIntersect.h @@ -0,0 +1,9 @@ +#ifndef TRITRIINTERSECT_H +#define TRITRIINTERSECT_H + +int tri_tri_intersect_with_isectline( + float V0[3],float V1[3],float V2[3], + float U0[3],float U1[3],float U2[3],int *coplanar, + float isectpt1[3],float isectpt2[3]); + +#endif |