Thanks for the interesting read!

Why is it that in your benchmark results, 1024 AABB tests are consistently faster than 32 AABB tests?

That's rather counterintuitive for me!

@Servant: I think the figures are per AABB (see the source code comments)

If you're using references, is it still a 'C implementation'?

Could you please choose better variable names than 'd' and 'r'? You then exacerbate the lack of clarity by aliasing d + r to 'd_p_r' and d - r to 'd_m_r' - not forgetting the 'd' field on the _Plane struct.

Would distanceFromCenterToPlaneNormal and distanceFromSizeToAbsolutePlaneNormal not be more explicative? Example code like this benefits hugely from such verbosity.

Thank you all for the comments.

@Servant: Bacterius is right. The numbers are "cycles per AABB". Culling a batch of 1024 AABBs in a single loop ends up being faster than 32 probably because the function overhead (e.g. calculating the abs of the plane normals) is minimal compared to the main loop.

Example: Assume that the initial loop which calculates the abs of the plane normals requires 200 cycles. Also, assume that each AABB requires 100 cycles. For a batch of 10 AABBs, the function would require 1200 cycles to complete. Or in other words, 120 cycles per AABB. If the batch had 1000 AABBs, the function would require 100200 cycles, or 100.2 cycles per AABB. Hope that makes sense.

@zdlr: Variable names were kept like that in order to match the examples in Fabian Giesen's article which I linked above. Also, the term 'reference implementation' doesn't mean that the code uses references. It means that, that specific snippet is used as a baseline for performance comparisons (if this was what you meant).

@hellraizer, thanks for replying. Re references, mean this:

const _Vector3f& aabbCenter = aabbList[iAABB].m_Center;

const _Vector3f& aabbSize = aabbList[iAABB].m_Extent;

Which creates constant references to the current AABB's center and extent vectors. References are a C++ construct, not C.

Your reference implementation can be optimized further.
You perform two dots per plane while it can be reduced to one.

Fabien Giesen covers this in excellent detail

On your 4 at a time SSE2 code, you're using movemask too often and switching to GPR registers to perform the bitwise arithmetic.
I suggest you keep everything in XMM registers (actually, avoid the switch of registers) and use cmpeqps familiy of instructions together with andps & orps instructions; then use movemask to get the final mask outside the loop

Another optimization is to analyze where to put the restrict/__restrict keyword to allow the compiler to optimize further.

@zdlr: Of course you are right. I forgot about those two. I'll edit the article to read "C++" instead of "C".

@Matias: Thank you very much for the tips. I'll try to find some time to do the changes you suggest and edit the article. I'll also try to make a little VS project to attach to the article at the same time.

About the two dot products. Do you mean method 5 on that article? In this case, the resulting code wouldn't be able to distinguish between fully inside and intersecting states which is a requirement for this article.

I know, it might sound bad trying to optimize something and add arbitrary restrictions which might affect performance, but I think having the ability to distinguish between those two cases might help in the case of an hierarchy. E.g. parent is completely inside, so there's no need to check any of its children.

Please correct me if I'm wrong.

Btw. a macro is often used to add the restrict keyword since every compiler has it's own way of declaring it (i.e. gcc vs msvc vs clang)

As for the two dot products, yes I'm talking about method 5. I wrongly saw in your code your routine was returning outside and inside and didn't see the "intersect". My apologies.

I wonder though if real life scenarios would just perform better using method 5 (brute force over smartness).

Also, it may be worth noting in the article that your routine is 100% re-entrant, meaning it's very easy to execute in parallel; except for your "static const unsigned int absPlaneMask[4]" part. That variable should be declared outside to make it re-entrant.

Being very strict, the snippet "_mm_load_ps((float*)&absPlaneMask[0])" is not standard compliant because it breaks the strict aliasing rule (casting an unsigned int* to a float*). It's a PITA coming from the intrinsics, since movaps can load integers just fine.

The "correct" (standard compliant) thing to do would be either to declare absPlaneMask as unsigned char* and then cast to float*

Another option is to rely on SSE2, and declare the variable as __m128 absPlaneMask( _mm_castsi128_ps( _mm_set1_epi32( 0x7FFFFFFF ) ) );

Absurd, but well, standard compliant (provided the compiler supports the SSE intrinsic extensions, which by definition aren't standard..., but you never know in which compiler your code will be used, and if it won't be ported to a different set of intrinsics)

Btw you shouldnt use _<Captial char> identifiers for C/C++ code as these are reserved as global names, this is one of the reasons why STL code looks so horrible. I am assuming here that this code should also compile for a C++ compiler.

17.4.3.2.1 Global names [lib.global.names]

Certain sets of names and function signatures are always reserved to the implementation:

• Each name that contains a double underscore (_ _) or begins with an underscore followed by an uppercase letter (2.11) is reserved to the implementation for any use.
• Each name that begins with an underscore is reserved to the implementation for use as a name in the global namespace.165

165) Such names are also reserved in namespace ::std (17.4.3.1).

Also defining a few macros for allignment and typedefing the intrinsics makes the code easier to read in the end.

@Matias Goldberg.

About the `static const unsigned int absPlaneMask[4]`. Could you explain how that stops it being reentrant? It is read-only and never modified.

Also, re-entrancy is not sufficient to make something thread-safe. It is possible to be re-entrant and yet not thread safe, and vice-versa.

@zdlr

The problem with static const inside a function scope is that initialization will happen the first time it is being used.

So, if two or more threads try to enter a function where there is a static variable (declared at that scope) for the first time, there's a race condition which can potentially leave the variable in an inconsistent state.

Once the function has been called once, this danger is no longer present and can be called from multiple threads safely.

In other words, if we look at the generated assembly, it's the variable is not really read-only. To ensure this doesn't happen we need to declare it somewhere else.

As for the re-entrancy, I didn't say it means it's thread safe, but rather makes it very easy to parallelize.

In this case, caller just needs to be sure it only gives the function access to a portion of the data, different to every thread, which is very straightforward to implement.

Perhaps the biggest thing to worry about is that the amount of data to be sent to each thread needs to be multiple of 4 (since the last snippet SSE is operating with 4 objects at a time)

Everyone, thanks for the constructive comments. I've updated the article based on some of your suggestions.

@Matias Goldberg: Unfortunately, changing _mm_load_ps((float*)&absPlaneMask[0]) to be standard compliant as you suggested, as well as adding the __restrict keyword, requires a rerun of the benchmarks, because otherwise the numbers won't be accurate. I'll keep a note and hope to be able to do it sometime soon.

One small note about your last comment. The 4 AABBs at a time SSE version has a loop at the end which should handle the case where the number of AABBs isn't a multiple of 4. The code isn't shown for clarity (it should actually be the same as the 1 AABB at a time version).

Thanks again for the great input.