Fixes and optimizations arount check_point_on_triangle functions

Author: sadko4u

CHANGELOG

@@ -3,6 +3,7 @@
 *******************************************************************************
 
 === 1.0.31 ===
+* Implementation of 3D mathematics using AVX instruction set.
 * Fixed improper AVX-512 optimization for lanczos kernel genration function.
 * Fixed Mid/Side transformation functions bug for 32-bit and 64-bit ARM processors. 
 * Added MacOS CI builds.

include/lsp-plug.in/dsp/common/3dmath.h

@@ -403,25 +403,29 @@ LSP_DSP_LIB_SYMBOL(float, check_triplet3d_vvn, const LSP_DSP_LIB_TYPE(vector3d_t
  */
 LSP_DSP_LIB_SYMBOL(float, check_triplet3d_vv, const LSP_DSP_LIB_TYPE(vector3d_t) *v);
 
-/** Analyze point location relative to the triangle of three points
+/**
+ * Analyze point location relative to the triangle of three points.
+ * All points should be considered being co-planar.
  *
  * @param t array of three triangle points
  * @param p point
- * @return value > 0 if point is candidate to be inside the triangle,
- * value < 0 if point is candidate to be outside the triangle,
- * value = 0 if point is on the edge of triangle
+ * @return value > 0 if point is inside of the triangle,
+ * value < 0 if point is outside of the triangle,
+ * value = 0 if point is on the edge of the triangle
  */
 LSP_DSP_LIB_SYMBOL(float, check_point3d_on_triangle_pvp, const LSP_DSP_LIB_TYPE(point3d_t) *t, const LSP_DSP_LIB_TYPE(point3d_t) *p);
 
-/** Analyze point location relative to the triangle of three points
+/**
+ * Analyze point location relative to the triangle of three points.
+ * All points should be considered being co-planar.
  *
  * @param p1 triangle point 1
  * @param p2 triangle point 2
  * @param p3 triangle point 3
  * @param p point
- * @return value > 0 if point is candidate to be inside the triangle,
- * value < 0 if point is candidate to be outside the triangle,
- * value = 0 if point is on the edge of triangle
+ * @return value > 0 if point is inside of the triangle,
+ * value < 0 if point is outside of the triangle,
+ * value = 0 if point is on the edge of the triangle
  */
 LSP_DSP_LIB_SYMBOL(float, check_point3d_on_triangle_p3p, const LSP_DSP_LIB_TYPE(point3d_t) *p1, const LSP_DSP_LIB_TYPE(point3d_t) *p2, const LSP_DSP_LIB_TYPE(point3d_t) *p3, const LSP_DSP_LIB_TYPE(point3d_t) *p);
 

include/private/dsp/arch/generic/3dmath.h

@@ -1112,16 +1112,7 @@ namespace lsp
  if (r[2] < 0.0f)
  return r[2];
 
- // Check 4
- r[2] = r[0]*r[1]*r[2];
- if (r[2] != 0.0f)
- return r[2];
-
- // Edge check: 3 scalar multiplications
- r[0] = v[0].dx * v[1].dx + v[0].dy * v[1].dy + v[0].dz * v[1].dz;
- r[1] = v[1].dx * v[2].dx + v[1].dy * v[2].dy + v[1].dz * v[2].dz;
- r[2] = v[2].dx * v[0].dx + v[2].dy * v[0].dy + v[2].dz * v[0].dz;
-
+ // The point may be located on the edge or on the vertex
  return r[0]*r[1]*r[2];
  }
 
@@ -1172,16 +1163,7 @@ namespace lsp
  if (r[2] < 0.0f)
  return r[2];
 
- // Check 4
- r[2] = r[0]*r[1]*r[2];
- if (r[2] != 0.0f)
- return r[2];
-
- // Edge check: 3 scalar multiplications
- r[0] = v[0].dx * v[1].dx + v[0].dy * v[1].dy + v[0].dz * v[1].dz;
- r[1] = v[1].dx * v[2].dx + v[1].dy * v[2].dy + v[1].dz * v[2].dz;
- r[2] = v[2].dx * v[0].dx + v[2].dy * v[0].dy + v[2].dz * v[0].dz;
-
+ // The point may be located on the edge or on the vertex
  return r[0]*r[1]*r[2];
  }
 

include/private/dsp/arch/x86/avx/3dmath.h

@@ -116,34 +116,34 @@
 
 /* 3x vector multiplication
  * Input:
- * x0 = vector1 [dx dy dz ? ]
- * x1 = vector2 [dx dy dz ? ]
- * x2 = vector3 [dx dy dz ? ]
+ * x0 = vector1 [dx dy dz 0 ]
+ * x1 = vector2 [dx dy dz 0 ]
+ * x2 = vector3 [dx dy dz 0 ]
  *
  * Output:
- * x0 = vector1 * vector2 [ vz vx vy ? ]
- * x1 = vector2 * vector3 [ vz vx vy ? ]
- * x2 = vector3 * vector1 [ vz vx vy ? ]
+ * x0 = vector1 * vector2 [ vz vx vy 0 ]
+ * x1 = vector2 * vector3 [ vz vx vy 0 ]
+ * x2 = vector3 * vector1 [ vz vx vy 0 ]
  */
 #define VECTOR_MUL3(x0, x1, x2, x3, x4, x5, x6, x7) \
- __ASM_EMIT("movaps %" x0 ", %" x3) /* xmm3 = dx0 dy0 dz0 dw0 */ \
- __ASM_EMIT("movaps %" x1 ", %" x4) /* xmm4 = dx1 dy1 dz1 dw1 */ \
- __ASM_EMIT("movaps %" x2 ", %" x5) /* xmm5 = dx2 dy2 dz2 dw2 */ \
- __ASM_EMIT("shufps $0xc9, %" x3 ", %" x3) /* xmm3 = dy0 dz0 dx0 dw0 */ \
- __ASM_EMIT("shufps $0xc9, %" x4 ", %" x4) /* xmm4 = dy1 dz1 dx1 dw1 */ \
- __ASM_EMIT("shufps $0xc9, %" x5 ", %" x5) /* xmm5 = dy2 dz2 dx2 dw2 */ \
- __ASM_EMIT("movaps %" x0 ", %" x6) /* xmm6 = dx0 dy0 dz0 dw0 */ \
- __ASM_EMIT("movaps %" x3 ", %" x7) /* xmm7 = dy0 dz0 dx0 dw0 */ \
+ __ASM_EMIT("movaps %" x0 ", %" x3) /* xmm3 = dx0 dy0 dz0 0 */ \
+ __ASM_EMIT("movaps %" x1 ", %" x4) /* xmm4 = dx1 dy1 dz1 0 */ \
+ __ASM_EMIT("movaps %" x2 ", %" x5) /* xmm5 = dx2 dy2 dz2 0 */ \
+ __ASM_EMIT("shufps $0xc9, %" x3 ", %" x3) /* xmm3 = dy0 dz0 dx0 0 */ \
+ __ASM_EMIT("shufps $0xc9, %" x4 ", %" x4) /* xmm4 = dy1 dz1 dx1 0 */ \
+ __ASM_EMIT("shufps $0xc9, %" x5 ", %" x5) /* xmm5 = dy2 dz2 dx2 0 */ \
+ __ASM_EMIT("movaps %" x0 ", %" x6) /* xmm6 = dx0 dy0 dz0 0 */ \
+ __ASM_EMIT("movaps %" x3 ", %" x7) /* xmm7 = dy0 dz0 dx0 0 */ \
  \
- __ASM_EMIT("mulps %" x4 ", %" x0) /* xmm0 = dx0*dy1 dy0*dz1 dz0*dx1 dw0*dw1 */ \
- __ASM_EMIT("mulps %" x1 ", %" x3) /* xmm3 = dy0*dx1 dz0*dy1 dx0*dz1 dw0*dw1 */ \
- __ASM_EMIT("mulps %" x5 ", %" x1) /* xmm1 = dx1*dy2 dy1*dz2 dz1*dx2 dw1*dw2 */ \
- __ASM_EMIT("mulps %" x2 ", %" x4) /* xmm4 = dy1*dx2 dz1*dy2 dx1*dz2 dw1*dw2 */ \
- __ASM_EMIT("mulps %" x7 ", %" x2) /* xmm2 = dx2*dy0 dy2*dz0 dz2*dx0 dw2*dw0 */ \
- __ASM_EMIT("mulps %" x6 ", %" x5) /* xmm5 = dy2*dx0 dz2*dy0 dx2*dz0 dw2*dw0 */ \
- __ASM_EMIT("subps %" x3 ", %" x0) /* xmm0 = nz0 nx0 ny0 nw0 */ \
- __ASM_EMIT("subps %" x4 ", %" x1) /* xmm1 = nz1 nx1 ny1 nw1 */ \
- __ASM_EMIT("subps %" x5 ", %" x2) /* xmm2 = nz2 nx2 ny2 nw2 */
+ __ASM_EMIT("mulps %" x4 ", %" x0) /* xmm0 = dx0*dy1 dy0*dz1 dz0*dx1 0 */ \
+ __ASM_EMIT("mulps %" x1 ", %" x3) /* xmm3 = dy0*dx1 dz0*dy1 dx0*dz1 0 */ \
+ __ASM_EMIT("mulps %" x5 ", %" x1) /* xmm1 = dx1*dy2 dy1*dz2 dz1*dx2 0 */ \
+ __ASM_EMIT("mulps %" x2 ", %" x4) /* xmm4 = dy1*dx2 dz1*dy2 dx1*dz2 0 */ \
+ __ASM_EMIT("mulps %" x7 ", %" x2) /* xmm2 = dx2*dy0 dy2*dz0 dz2*dx0 0 */ \
+ __ASM_EMIT("mulps %" x6 ", %" x5) /* xmm5 = dy2*dx0 dz2*dy0 dx2*dz0 0 */ \
+ __ASM_EMIT("subps %" x3 ", %" x0) /* xmm0 = nz0 nx0 ny0 0 */ \
+ __ASM_EMIT("subps %" x4 ", %" x1) /* xmm1 = nz1 nx1 ny1 0 */ \
+ __ASM_EMIT("subps %" x5 ", %" x2) /* xmm2 = nz2 nx2 ny2 0 */
 
 /* 1x vector multiplication
  * Input:
@@ -1331,49 +1331,39 @@ namespace lsp
  ARCH_X86_ASM
  (
  /* Load vectors */
- __ASM_EMIT("movups (%[p]), %[x3]") /* xmm3 = px py pz pw */
- __ASM_EMIT("movups 0x00(%[pv]), %[x0]") /* xmm0 = x0 y0 z0 w0 */
- __ASM_EMIT("movups 0x10(%[pv]), %[x1]") /* xmm1 = x1 y1 z1 w1 */
- __ASM_EMIT("movups 0x20(%[pv]), %[x2]") /* xmm2 = x2 y2 z2 w2 */
- __ASM_EMIT("subps %[x3], %[x0]") /* xmm0 = dx0 dy0 dz0 dw0 */
- __ASM_EMIT("subps %[x3], %[x1]") /* xmm1 = dx1 dy1 dz1 dw1 */
- __ASM_EMIT("subps %[x3], %[x2]") /* xmm2 = dx2 dy2 dz2 dw2 */
+ __ASM_EMIT("movups (%[p]), %[x3]") /* xmm3 = px py pz 1 */
+ __ASM_EMIT("movups 0x00(%[pv]), %[x0]") /* xmm0 = x0 y0 z0 1 */
+ __ASM_EMIT("movups 0x10(%[pv]), %[x1]") /* xmm1 = x1 y1 z1 1 */
+ __ASM_EMIT("movups 0x20(%[pv]), %[x2]") /* xmm2 = x2 y2 z2 1 */
+ __ASM_EMIT("subps %[x3], %[x0]") /* xmm0 = v0 = dx0 dy0 dz0 0 */
+ __ASM_EMIT("subps %[x3], %[x1]") /* xmm1 = v1 = dx1 dy1 dz1 0 */
+ __ASM_EMIT("subps %[x3], %[x2]") /* xmm2 = v2 = dx2 dy2 dz2 0 */
  /* 3x vector multiplications */
  VECTOR_MUL3("[x0]", "[x1]", "[x2]", "[x3]", "[x4]", "[x5]", "[x6]", "[x7]")
  /* 3x scalar multiplications */
+ /* xmm0 = m0 = v0 cross v1 */
+ /* xmm1 = m1 = v1 cross v2 */
+ /* xmm2 = m2 = v2 cross v0 */
  SCALAR_MUL3("[x0]", "[x1]", "[x2]", "[x3]")
  /* Compare with zeros */
- __ASM_EMIT("xorps %[x4], %[x4]")
- __ASM_EMIT("ucomiss %[x4], %[x0]")
- __ASM_EMIT("jb 110f")
- __ASM_EMIT("ucomiss %[x4], %[x1]")
- __ASM_EMIT("jb 109f")
- __ASM_EMIT("ucomiss %[x4], %[x2]")
- __ASM_EMIT("jb 108f")
- __ASM_EMIT("mulss %[x1], %[x0]")
- __ASM_EMIT("mulss %[x2], %[x0]")
- __ASM_EMIT("ucomiss %[x4], %[x0]")
- __ASM_EMIT("jne 110f")
- /* There is somewhere zero, need additional check */
- /* Load vectors */
- __ASM_EMIT("movups (%[p]), %[x3]") /* xmm3 = px py pz pw */
- __ASM_EMIT("movups 0x00(%[pv]), %[x0]") /* xmm0 = x0 y0 z0 w0 */
- __ASM_EMIT("movups 0x10(%[pv]), %[x1]") /* xmm1 = x1 y1 z1 w1 */
- __ASM_EMIT("movups 0x20(%[pv]), %[x2]") /* xmm2 = x2 y2 z2 w2 */
- __ASM_EMIT("subps %[x3], %[x0]") /* xmm0 = dx0 dy0 dz0 dw0 */
- __ASM_EMIT("subps %[x3], %[x1]") /* xmm1 = dx1 dy1 dz1 dw1 */
- __ASM_EMIT("subps %[x3], %[x2]") /* xmm2 = dx2 dy2 dz2 dw2 */
- /* Do 3x scalar multiplications */
- SCALAR_MUL3("[x0]", "[x1]", "[x2]", "[x3]")
- __ASM_EMIT("mulss %[x1], %[x0]")
- __ASM_EMIT("mulss %[x2], %[x0]")
- __ASM_EMIT("jmp 110f")
-
- __ASM_EMIT("108:")
- __ASM_EMIT("movss %[x2], %[x0]")
+ /* xmm0 = r0 = m0 dot m1 */
+ /* xmm1 = r1 = m1 dot m2 */
+ /* xmm2 = r2 = m2 dot m0 */
+ __ASM_EMIT("xorps %[x4], %[x4]") /* xmm4 = 0 */
+ __ASM_EMIT("ucomiss %[x4], %[x0]") /* r0 <=> 0 */
+ __ASM_EMIT("jb 110f") /* r0 < 0 */
+ __ASM_EMIT("mulss %[x1], %[x0]") /* xmm0 = r0 * r1 */
+ __ASM_EMIT("ucomiss %[x4], %[x1]") /* r1 <=> 0 */
+ __ASM_EMIT("jb 109f") /* r1 < 0 */
+ __ASM_EMIT("mulss %[x2], %[x0]") /* xmm0 = r0 * r1 * r2 */
+ __ASM_EMIT("ucomiss %[x4], %[x2]") /* r2 <=> 0 */
+ __ASM_EMIT("jae 110f") /* r2 >= 0 */
+ /* Fail cases */
+ __ASM_EMIT("movaps %[x2], %[x0]")
  __ASM_EMIT("jmp 110f")
  __ASM_EMIT("109:")
- __ASM_EMIT("movss %[x1], %[x0]")
+ __ASM_EMIT("movaps %[x1], %[x0]")
+ /* End */
  __ASM_EMIT("110:")
  : [x0] "=&x" (x0), [x1] "=&x" (x1), [x2] "=&x" (x2), [x3] "=&x" (x3),
  [x4] "=&x" (x4), [x5] "=&x" (x5), [x6] "=&x" (x6), [x7] "=&x" (x7)
@@ -1390,49 +1380,39 @@ namespace lsp
  ARCH_X86_ASM
  (
  /* Load vectors */
- __ASM_EMIT("movups (%[p]), %[x3]") /* xmm3 = px py pz pw */
- __ASM_EMIT("movups (%[p1]), %[x0]") /* xmm0 = x0 y0 z0 w0 */
- __ASM_EMIT("movups (%[p2]), %[x1]") /* xmm1 = x1 y1 z1 w1 */
- __ASM_EMIT("movups (%[p3]), %[x2]") /* xmm2 = x2 y2 z2 w2 */
- __ASM_EMIT("subps %[x3], %[x0]") /* xmm0 = dx0 dy0 dz0 dw0 */
- __ASM_EMIT("subps %[x3], %[x1]") /* xmm1 = dx1 dy1 dz1 dw1 */
- __ASM_EMIT("subps %[x3], %[x2]") /* xmm2 = dx2 dy2 dz2 dw2 */
+ __ASM_EMIT("movups (%[p]), %[x3]") /* xmm3 = px py pz 1 */
+ __ASM_EMIT("movups (%[p1]), %[x0]") /* xmm0 = x0 y0 z0 1 */
+ __ASM_EMIT("movups (%[p2]), %[x1]") /* xmm1 = x1 y1 z1 1 */
+ __ASM_EMIT("movups (%[p3]), %[x2]") /* xmm2 = x2 y2 z2 1 */
+ __ASM_EMIT("subps %[x3], %[x0]") /* xmm0 = v0 = dx0 dy0 dz0 0 */
+ __ASM_EMIT("subps %[x3], %[x1]") /* xmm1 = v1 = dx1 dy1 dz1 0 */
+ __ASM_EMIT("subps %[x3], %[x2]") /* xmm2 = v2 = dx2 dy2 dz2 0 */
  /* 3x vector multiplications */
  VECTOR_MUL3("[x0]", "[x1]", "[x2]", "[x3]", "[x4]", "[x5]", "[x6]", "[x7]")
  /* 3x scalar multiplications */
+ /* xmm0 = m0 = v0 cross v1 */
+ /* xmm1 = m1 = v1 cross v2 */
+ /* xmm2 = m2 = v2 cross v0 */
  SCALAR_MUL3("[x0]", "[x1]", "[x2]", "[x3]")
  /* Compare with zeros */
- __ASM_EMIT("xorps %[x4], %[x4]")
- __ASM_EMIT("ucomiss %[x4], %[x0]")
- __ASM_EMIT("jb 110f")
- __ASM_EMIT("ucomiss %[x4], %[x1]")
- __ASM_EMIT("jb 109f")
- __ASM_EMIT("ucomiss %[x4], %[x2]")
- __ASM_EMIT("jb 108f")
- __ASM_EMIT("mulss %[x1], %[x0]")
- __ASM_EMIT("mulss %[x2], %[x0]")
- __ASM_EMIT("ucomiss %[x4], %[x0]")
- __ASM_EMIT("jne 110f")
- /* There is somewhere zero, need additional check */
- /* Load vectors */
- __ASM_EMIT("movups (%[p]), %[x3]") /* xmm3 = px py pz pw */
- __ASM_EMIT("movups (%[p1]), %[x0]") /* xmm0 = x0 y0 z0 w0 */
- __ASM_EMIT("movups (%[p2]), %[x1]") /* xmm1 = x1 y1 z1 w1 */
- __ASM_EMIT("movups (%[p3]), %[x2]") /* xmm2 = x2 y2 z2 w2 */
- __ASM_EMIT("subps %[x3], %[x0]") /* xmm0 = dx0 dy0 dz0 dw0 */
- __ASM_EMIT("subps %[x3], %[x1]") /* xmm1 = dx1 dy1 dz1 dw1 */
- __ASM_EMIT("subps %[x3], %[x2]") /* xmm2 = dx2 dy2 dz2 dw2 */
- /* Do 3x scalar multiplications */
- SCALAR_MUL3("[x0]", "[x1]", "[x2]", "[x3]")
- __ASM_EMIT("mulss %[x1], %[x0]")
- __ASM_EMIT("mulss %[x2], %[x0]")
- __ASM_EMIT("jmp 110f")
-
- __ASM_EMIT("108:")
- __ASM_EMIT("movss %[x2], %[x0]")
+ /* xmm0 = r0 = m0 dot m1 */
+ /* xmm1 = r1 = m1 dot m2 */
+ /* xmm2 = r2 = m2 dot m0 */
+ __ASM_EMIT("xorps %[x4], %[x4]") /* xmm4 = 0 */
+ __ASM_EMIT("ucomiss %[x4], %[x0]") /* r0 <=> 0 */
+ __ASM_EMIT("jb 110f") /* r0 < 0 */
+ __ASM_EMIT("mulss %[x1], %[x0]") /* xmm0 = r0 * r1 */
+ __ASM_EMIT("ucomiss %[x4], %[x1]") /* r1 <=> 0 */
+ __ASM_EMIT("jb 109f") /* r1 < 0 */
+ __ASM_EMIT("mulss %[x2], %[x0]") /* xmm0 = r0 * r1 * r2 */
+ __ASM_EMIT("ucomiss %[x4], %[x2]") /* r2 <=> 0 */
+ __ASM_EMIT("jae 110f") /* r2 >= 0 */
+ /* Fail cases */
+ __ASM_EMIT("movaps %[x2], %[x0]")
  __ASM_EMIT("jmp 110f")
  __ASM_EMIT("109:")
- __ASM_EMIT("movss %[x1], %[x0]")
+ __ASM_EMIT("movaps %[x1], %[x0]")
+ /* End */
  __ASM_EMIT("110:")
  : [x0] "=&x" (x0), [x1] "=&x" (x1), [x2] "=&x" (x2), [x3] "=&x" (x3),
  [x4] "=&x" (x4), [x5] "=&x" (x5), [x6] "=&x" (x6), [x7] "=&x" (x7)

include/private/dsp/arch/x86/sse/3dmath.h

@@ -488,6 +488,9 @@
  EXPORT1(check_triplet3d_vvn);
  EXPORT1(check_triplet3d_vv);
 
+ EXPORT1(check_point3d_on_triangle_p3p);
+ EXPORT1(check_point3d_on_triangle_pvp);
+
  // FMA3 support?
  if (f->features & CPU_OPTION_FMA3)
  {
@@ -609,6 +612,9 @@
  CEXPORT2(favx, check_triplet3d_v2n, check_triplet3d_v2n_fma3);
  CEXPORT2(favx, check_triplet3d_vvn, check_triplet3d_vvn_fma3);
  CEXPORT2(favx, check_triplet3d_vv, check_triplet3d_vv_fma3);
+
+ CEXPORT2(favx, check_point3d_on_triangle_p3p, check_point3d_on_triangle_p3p_fma3);
+ CEXPORT2(favx, check_point3d_on_triangle_pvp, check_point3d_on_triangle_pvp_fma3);
  }
  }