Jump to content

  • Log In with Google      Sign In   
  • Create Account

Dave Eberly

Member Since 22 Aug 2004
Offline Last Active Oct 15 2012 01:13 AM

#4977502 Figuring out ellipsoid-ellipsoid collision detection

Posted by on 06 September 2012 - 11:50 PM

My advice is to skip ellipsoids. Use a bounding box or a k-DOP (or some convex polyhedron with a small number of faces), then use separating axis tests. The coding is just a lot easier, and the numerical robustness issues in determining ellipsoid-ellipsoid intersection can be avoided.

That said, there is an implementation of the ellipsoid-ellipsoid intersection at my web site.

Regarding capsules, better choice than ellipsoids because the intersection/separation tests are simpler to implement. For capsule-capsule sweep, a simple implementation uses bisection over the desired time interval. At each time you apply a (static) capsule-capsule overlap test, which reduces to a computation of distance between line segments and a comparson involving this distance and capsule radii. You can avoid the iteration--I have pseudocode for this in my Game Physics 2nd edition (section 6.3.2).

Regarding bounded cylinders, the game physics book and a document at my web site show the complexity of intersection testing via separating axes. Turns out that it is simpler to do cylinder-cylinder intersection with infinite cylinders, then clip the intersection set based on the finite cylinder heights. Not an exact test (result is not a closed form solution), but effective. I have a document about this at my web site and a sample application (for the infinite cylinder-cylinder test).

#4968344 Hieroglyph 3 Rendering engine Question

Posted by on 10 August 2012 - 11:46 PM

Is the Hieroglyph 3 well written?
I am trying to learn the right way to make dynamic classes for a engine like rendering system using direct3d 11.

Is there another good open source engine using direct3d?

I want to know how the professionals make the thinks.

The code is good quality. More importantly, purchase the book that goes with it: Practical Rendering & Computation with Direct3D 11, by J. Zink, M. Pettineo, and J. Hoxley. I have requested reading this for the engineers on my real-time graphics team.

#4965099 SSE vector normalization

Posted by on 31 July 2012 - 11:44 PM

inline const CVector3SSE& CVector3SSE::Normalize()
  static const __m128 almostZero = _mm_set1_ps(1e-5f);
  __m128 dp = _mm_dp_ps(m_fValsSSE, m_fValsSSE, 0x7F);
  const __m128 cmp = _mm_gt_ps(dp, almostZero);
  dp = _mm_rsqrt_ps(dp);
  m_fValsSSE = _mm_mul_ps(m_fValsSSE, _mm_and_ps(dp, cmp));
  return *this;

Although yours is the standard way folks do the normalization, for large components the dot product overflows. If you need something that is robust for all finite floating-point inputs,
inline __m128 MaximumAbsoluteComponent (__m128 const v)
	__m128 SIGN = _mm_set1_ps(0x80000000u);
	__m128 vAbs = _mm_andnot_ps(SIGN, v);
	__m128 max0 = _mm_shuffle_ps(vAbs, vAbs, _MM_SHUFFLE(0,0,0,0));
	__m128 max1 = _mm_shuffle_ps(vAbs, vAbs, _MM_SHUFFLE(1,1,1,1));
	__m128 max2 = _mm_shuffle_ps(vAbs, vAbs, _MM_SHUFFLE(2,2,2,2));
	__m128 max3 = _mm_shuffle_ps(vAbs, vAbs, _MM_SHUFFLE(3,3,3,3));
	max0 = _mm_max_ps(max0, max1);
	max2 = _mm_max_ps(max2, max3);
	max0 = _mm_max_ps(max0, max2);
	return max0;

inline __m128 Normalize (__m128 const v)
	// Compute the maximum absolute value component.
	__m128 maxComponent = MaximumAbsoluteComponent(v);

	// Divide by the maximum absolute component.  This is potentially a divide by zero.
	__m128 normalized = _mm_div_ps(v, maxComponent);

	// Set to zero when the original length is zero.
	__m128 zero = _mm_setzero_ps();
	__m128 mask = _mm_cmpneq_ps(zero, maxComponent);
	normalized = _mm_and_ps(mask, normalized);

	// (sqrLength, sqrLength, sqrLength, sqrLength)
	__m128 sqrLength = _mm_dp_ps(normalized, normalized, 0x7F);

	// (length, length, length, length)
	__m128 length = _mm_sqrt_ps(sqrLength);

	// Divide by the length to normalize.  This is potentially a divide by zero.
	normalized = _mm_div_ps(normalized, length);

	// Set to zero when the original length is zero or infinity.  In the latter case, this is considered to be an unexpected condition.
	normalized = _mm_and_ps(mask, normalized);
	return normalized;

#4930121 Trying to find a minimal set of polygons that intersect a rectangle.

Posted by on 10 April 2012 - 11:27 PM

I have an ordered set of polygons and I am intersecting them with a rectangle, how do I tell when the rectangle is filled so I can stop checking? This is in 2D.

I am trying to get the minimal set of polygons with cover the rectangle. The order of the set is the priority of use and the polygons are simple (don't have holes or self intersecting). So if my set is {A,B,C,D} and A B D intersect the rectangle but A contains B than I want a the set that contains just {A, D}.

Tried googling but wasn't sure what to google really, I was thinking it would fall under polygon coverage but that didn't seem to produce results of what I was looking for.


I think this is a hard problem theoretically. For a practical solution, have you thought about rasterizing the rectangle and polygons to a high-resolution grid? Rasterize the rectangle first. Rasterizer your polygons one at a time, keeping track of which rectangle pixels are written to (sort of like a stencil buffer) and which polygons have been rasterized to previously unwritten pixels. Once you have rasterized to all pixels, the process terminates. (You have to deal with not writing all rectangle pixels after all polygons have been processed.)

#4919702 SLMATH library and SSE optimisation problem.

Posted by on 06 March 2012 - 01:51 AM

Thanks. It does, yes. It's a giant pain, so I think I'll just switch off SSE in the project settings and live without the performance boost.

std::vector should be able to support alignment through custom allocators. However, if you are using MSVS 2010, the dinkumware STL they use has a bug in that the std::vector resize does not do the right thing (fixed in MSVS 2011). For MSVS 2010, you'll have to roll your own std::vector (maybe copy what dinkumware does and "fix" the resize).

#4915790 ID3D11ShaderReflection question

Posted by on 23 February 2012 - 12:16 AM

I found another post of yours that mentions the assignment of each element of a float[] to a single register. That was helpful. When I query for the array member, I saw Rows=1, Columns=1, and Elements=5, which seemed strange (a 1x1 array with 5 elements?). Because the query does not tell me I have an "array", I suppose that I can infer it from Rows==1 and Columns==1 and Elements > 1. The other information that seemed strange: cbuffer Whatever { struct Something {...}; Something A; Something B[2]}. The number of bytes reported for A is different from the number of bytes reported for B. Of course, for B the number of bytes is for both array items, but it was not twice the size of A. And it appears I'd have to assume that both array items of B use the same number of bytes. Reverse engineering (by compiling some shaders that access the .x components of various struct members) made it clear that there was some set of rules the compiler was using. So I think that knowing how arrays are mapped to registers and knowing to trap the rows/columns/elements case I mentioned, I can infer the packing of the struct (which I want so I know how to typecast the mapped cbuffer data properly). Thanks.

#4860227 Sweep Line Algorithm

Posted by on 10 September 2011 - 11:03 PM

What you suggest works for the initial sort of the boxes. However, when you start moving the boxes and updating each frame, a sort only on x is not sufficient. Consider box0 with 0 <= x <= 2 and 0 <= y <= 2. Consider box1 with 1 <= x <= 3 and 1 <= y <= 3. In the initial sort on x, you find that x-intervals [0,2] and [1,3] overlap, so then you test for box-box overlap and find the boxes intersect. Suppose box1 is moved to 1 <= x <= 3 and 3 <= y <= 4; that is, the box is moved in the y-direction. The two boxes are no longer overlapping, so you need to remove the pair from the overlap set. If you only sort on x, you would never enter the insertion-sort swap code because the x-endpoints have not changed at all (the x-list is still correctly ordered).

#4847025 OBB Computation

Posted by on 09 August 2011 - 10:36 PM

If you want as good a fit as possible, the minimum-volume OBB is a decent choice. As mentioned, it is slow, so probably better for off-line computations rather than during the real-time application.

Alternatively, you can use the covariance approach. If you have two (nearly) equal eigenvalues (eigenspace of dimension 2) and a third eigenvalue that is different from them (eigenspace of dimension 1), you can project out the eigenvector from the third eigenvalue. You now have a model in the 2-dimensional eigenspace. Now you can use the method of rotating calipers to find the minimum-area rectangle containing the (projected) points. Effectively, this allows you to search the 2-dimensional eigenspace for a pair of orthogonal eigenvectors that give you a good fit. I suspect this approach is faster than always trying to compute the minimum-volume OBB.

#4846520 OBB Computation

Posted by on 08 August 2011 - 10:09 PM

The problem with a square is that the covariance matrix has an eigenvalue of multiplicity 2; that is, the dimension of the eigenspace is 2, so there are infinitely many choices for a pair of orthogonal eigenvectors. A numerical eigensolver is free to compute any of these pairs, and you cannot expect that the pair it chooses leads to the "best fit" you are expecting.. Generally, if there is a symmetry in the input data, you will probably have eigenspaces of dimension 2 or 3, and the OBB might not fit well for the eigenvectors returned by the numerical eigensolver.

In your square example in 3D, the center does not look right. The center should be the average of the vertices, which is (0,0,1), not (-0.5,-0.5,1). So your implementation might have other problems. Compute the center first. When building the covariance matrix, do not forget to subtract the center from each vertex before computing the sum of squared values.

#4841423 Closest point on ray to AABB?

Posted by on 27 July 2011 - 09:21 PM

See this page , section "Distance between lines, rays, or segments and boxes (3D). "

#4837794 Math.Sin replaced with lookup table

Posted by on 19 July 2011 - 10:48 PM

I am working currently on a core 2 duo in my macbook, but I am targeting the xbox 360 with XNA. I suspect that I am frequently accessing RAM because I am filling 9 buffers 48000 floats long. I am not too sure what the profiler is telling me because I am very inexperienced with it. It is the one that comes with visual studio 2010 professional.

I assume you are using intrinsics to get the SIMD versions of the XNA math functions? The XNA sine function uses a Taylor series with 12 terms. Using a different algorithm (Chebyshev Equioscillation Theorem applied to sin(x)/x), you can double the speed (with a polynomial of 6 terms) at half the maximum error. From my math libraries. It is easy to convert this to SIMD, whether Intel SSE or Xbox360.

template <typename Real>
Real Math<Real>::FastSin1 (Real angle)
    Real angleSqr = angle*angle;
    Real result = -(Real)2.39e-08;
    result *= angleSqr;
    result += (Real)2.7526e-06;
    result *= angleSqr;
    result -= (Real)1.98409e-04;
    result *= angleSqr;
    result += (Real)8.3333315e-03;
    result *= angleSqr;
    result -= (Real)1.666666664e-01;
    result *= angleSqr;
    result += (Real)1.0;
    result *= angle;
    return result;

#4834630 3D FingerPrint - need to extract ridges and valleys? Estimation Curvatures

Posted by on 12 July 2011 - 11:18 PM

You posted a similar request to this forum on June 22. I posted a link to a PDF about ridges/valleys and a link to a sample application. Did you try these ideas? If so, what did not work for you? You posted code, but folks here generally are reluctant to critique large code blocks, prefering instead to discuss your ideas.

You wrote: "I am totally frustated !! Tried all done all situation but nothing is working for me. I have to complete this within one month to defend my thesis. If I am not able to do so I will lose my job that I want to join after that "

The topic area is a difficult one. If you have made no progress at this time, one more month is not going to make a difference. Writing a thesis is one thing, defending it is another, not everyone successfully defends, and whether or not you obtain a job is not relevant to folks helping you here (and is not an incentive to make them want to help). Perhaps now is a good time to meet with your advisor and discuss your situation...

#4828953 Minimum bounding sphere of a Frustum

Posted by on 28 June 2011 - 11:06 PM


void ComputeMinimumVolumeSphereOfFrustum (const float eye[3],

    const float direction[3], const float dmin, const float dmax,

    const float rmax, const float umax, float center[3], float& radius)


    float invDMin = 1.0f/dmin;

    float u = umax*invDMin;

    float r = rmax*invDMin;

    float u2pr2 = u*u + r*r;

    float s = 0.5f*(dmin + dmax)*(1.0f + u2pr2);

    if (s >= dmax)


        s = dmax;

        radius = dmax*sqrt(u2pr2);




        float diff = 1.0f - s/dmax;

        radius = dmax*sqrt(diff*diff + u2pr2);


    for (int i = 0; i < 3; ++i)


        center[i] = eye[i] + s*direction[i];




template <int n, typename Real>

void ComputeMiniballFrustum (const Real eye[n], const Real direction[n],

    const Real min0, const Real max[n], Real center[n], Real& radius)


    Real invMin0 = ((Real)1)/min0;

    float ratio[n], sqrLength = (Real)0;

    for (int i = 1; i < n; ++i)


        ratio[i] = max[i]*invMin0;

        sqrLength += ratio[i]*ratio[i];


    Real s = ((Real)0.5)*(min0 + max[0])*((Real)1 + sqrLength);

    if (s >= max[0])


        s = max[0];

        radius = max[0]*sqrt(sqrLength);




        Real diff = (Real)1 - s/max[0];

        radius = max[0]*sqrt(diff*diff + sqrLength);


    for (int i = 0; i < n; ++i)


        center[i] = eye[i] + s*direction[i];




int main ()


    float eye[3] = { 0.0f, 0.0f, 0.0f };

    float direction[3] = { 0.0f, 0.0f, 1.0f };

    float dmin = 10.0f, dmax = 1000.0f, rmax = 2.0f, umax = 1.0f;

    float center0[3], radius0;

    ComputeMinimumVolumeSphereOfFrustum(eye, direction, dmin,

        dmax, rmax, umax, center0, radius0);

    float min0 = dmin, max[3] = { dmax, rmax, umax };

    float center1[3], radius1;

    ComputeMiniballFrustum<3,float>(eye, direction, min0, max, center1,


    return 0;



#4827954 Minimum bounding sphere of a Frustum

Posted by on 26 June 2011 - 11:43 AM

Let E be the eyepoint. Let the view direction be D, the up direction be U, and the right direction be R. D, U, and R are unit length and mutually perpendicular. The frustum values are dmin (near), dmax (far), rmax (right), -rmax (left), umax (up), and -umax (down). The sphere center is C = E + s*D for some s. A vertex of the near plane is P0 = E + dmin*D + umax*U + rmax*R. A vertex of the far plane is P1 = E + dmax*D + (dmax*umax/dmin)*U + (dmax*rmax/dmin)*R. The squared distances from P0 to C and P1 to C must be equal, so (s-dmin)^2 + rmax^2 + umax^2 = (s - dmax)^2 + (dmax*rmax/dmin)^2 + (dmax*umax/dmin)^2. This leads to s = ((dmin + dmax)/2)*(1 + (rmax/dmin)^2 + (umax/dmin)^2).

This is actually the circumscribed sphere, which I believe is the minimum-volume sphere as long as center C is inside the frustum, which means as long as s <= dmax. When s > dmax, the minimum-volume sphere is the one whose center is E+dmax*D. The radius is distance from C to a vertex on the far plane. In the case s < dmax, the radius is the distance from C to any of the 8 vertices, but you can just substitute the value of s in either of the squared distances I mentioned and then take the square root.

#4806750 Closest point on cubic bézier curve

Posted by on 04 May 2011 - 10:51 PM

I'm trying to find the point on a cubic bézier curve closest to another external point P, in 2D space. I've read the Graphics Gems explanation on the method and other stuff but the math's way over my head. Specifically I don't really understand how finding the roots of the.. curve's polynomial (am I making sense?) gets you closer to the solution. Would somebody be able to explain this in a 14 year old's terms? (:

Generally, let the curve be parameterized by (x(t),y(t)) and let P = (p0,p1) be the external point. The vector from the closest point on the curve to P must be normal to the curve, which means that vector is perpendicular to a tangent of the curve. Such a tangent is the derivative (x'(t),y'(t)). The vector from a curve point to P is (x(t)-p0,y(t)-p1). To be perpendicular, the dot product is zero: 0 = x'(t)*(x(t)-p0) + y'(t)*(y(t)-p1) = f(t). When (x(t),y(t)) has polynomial components, f(t) is a polynomial. So the problem reduces to finding roots of f(t). In your case, x(t) and y(t) are degree-3 polynomials, x'(t) and y'(t) are degree-2 polynomials, so f(t) is a degree-5 polynomial. There are no closed-form equations for roots of degree-5, so you must use a numerical method for root finding.