culling.c

// Brute-force call using world-space bounding spheres.
for (buc_uint32_t frustum = 0; frustum < n_frusta; ++frustum) {
  __m128 p_x0x1x2x3; __m128 p_y0y1y2y3; __m128 p_z0z1z2z3; __m128 p_d0d1d2d3;
  __m128 p_x4x5x4x5; __m128 p_y4y5y4y5; __m128 p_z4z5z4z5; __m128 p_d4d5d4d5;

  // Transpose frustum planes so we can test against them in parallel.
  {
    const __m128 x0y0z0d0 = _mm_load_ps((const float *)&frusta[frustum].planes[0]);
    const __m128 x1y1z1d1 = _mm_load_ps((const float *)&frusta[frustum].planes[1]);
    const __m128 x2y2z2d2 = _mm_load_ps((const float *)&frusta[frustum].planes[2]);
    const __m128 x3y3z3d3 = _mm_load_ps((const float *)&frusta[frustum].planes[3]);
    const __m128 x4y4z4d4 = _mm_load_ps((const float *)&frusta[frustum].planes[4]);
    const __m128 x5y5z5d5 = _mm_load_ps((const float *)&frusta[frustum].planes[5]);

    const __m128 x0x1y0y1 = _mm_unpackhi_ps(x1y1z1d1, x0y0z0d0);
    const __m128 x2x3y2y3 = _mm_unpackhi_ps(x3y3z3d3, x2y2z2d2);
    const __m128 z0z1d0d1 = _mm_unpacklo_ps(x1y1z1d1, x0y0z0d0);
    const __m128 z2z3d2d3 = _mm_unpacklo_ps(x3y3z3d3, x2y2z2d2);

    p_x0x1x2x3 = _mm_movehl_ps(x0x1y0y1, x2x3y2y3);
    p_y0y1y2y3 = _mm_movelh_ps(x0x1y0y1, x2x3y2y3);
    p_z0z1z2z3 = _mm_movehl_ps(z0z1d0d1, z2z3d2d3);
    p_d0d1d2d3 = _mm_movehl_ps(z0z1d0d1, z2z3d2d3);

    const __m128 x4x5y4y5 = _mm_unpackhi_ps(x5y5z5d5, x4y4z4d4);
    const __m128 z4z5d4d5 = _mm_unpacklo_ps(x5y5z5d5, x4y4z4d4);

    p_x4x5x4x5 = _mm_movehl_ps(x4x5y4y5, x4x5y4y5);
    p_y4y5y4y5 = _mm_movelh_ps(x4x5y4y5, x4x5y4y5);
    p_z4z5z4z5 = _mm_movehl_ps(z4z5d4d5, z4z5d4d5);
    p_d4d5d4d5 = _mm_movelh_ps(z4z5d4d5, z4z5d4d5);
  }

  // Build look up table mapping from test results into flag.
  const buc_uint32_t mask_to_flag[16] =
    { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, (1 << frustum) };

  // Test every entity against this frustum.
  for (buc_uint32_t i = 0; i < manager->n_entities; ++i) {
    const buc_entity_t *entity = &manager->entities[i];

    // Derive world-space bounding sphere.
    __m128 s_xxxx = _mm_set1_ps(entity->position.x);
    __m128 s_yyyy = _mm_set1_ps(entity->position.y);
    __m128 s_zzzz = _mm_set1_ps(entity->position.z);
    __m128 s_rrrr = _mm_set1_ps(entity->radius * ws_magnitude(entity->scale));

    // Compute dot products against planes.
    __m128 dot_0123 = _mm_mul_ps(s_xxxx, p_x0x1x2x3);
           dot_0123 = _mm_add_ps(_mm_mul_ps(s_yyyy, p_y0y1y2y3), dot_0123);
           dot_0123 = _mm_add_ps(_mm_mul_ps(s_zzzz, p_z0z1z2z3), dot_0123);
    __m128 dot_4545 = _mm_mul_ps(s_xxxx, p_x4x5x4x5);
           dot_4545 = _mm_add_ps(_mm_mul_ps(s_yyyy, p_y4y5y4y5), dot_4545);
           dot_4545 = _mm_add_ps(_mm_mul_ps(s_zzzz, p_z4z5z4z5), dot_4545);

    // Compute distances to planes.
    __m128 dist_0123 = _mm_sub_ps(dot_0123, p_d0d1d2d3);
    __m128 dist_4545 = _mm_sub_ps(dot_4545, p_d4d5d4d5);

    // Check if at least partially in positive half-space of each plane.
    __m128 results = _mm_cmpge_ps(dist_0123, s_rrrr);
           results = _mm_and_ps(results, _mm_cmpge_ps(dist_4545, s_rrrr));

    // The sphere is at least partially inside the frustum if and only if
    // it is at least partially in the half-space of every plane.
    const buc_uint32_t flag = mask_to_flag[_mm_movemask_ps(results)];

    // Set flag if at least partially in this frustum.
    visibility[i] |= flag;
  }
}