docs: @Pure
  	public static boolean containsTrianglePoint(
  			double ax, double ay, double az,
  			double bx, double by, double bz,
  			double cx, double cy, double cz,
  			double px, double py, double pz,
  			boolean forceCoplanar, double epsilon) {
  
  		//
  		// Compute vectors        
  		//
  		// v0 = C - A
  		double v0x = cx - ax;
  		double v0y = cy - ay;
  		double v0z = cz - az;
  		// v1 = B - A
  		double v1x = bx - ax;
  		double v1y = by - ay;
  		double v1z = bz - az;
  		// v2 = P - A
  		double v2x = px - ax;
  		double v2y = py - ay;
  		double v2z = pz - az;
  
  		//
  		// Compute dot products
  		//
  		// dot01 = dot(v0, v0)
  		double dot00 = FunctionalVector3D.dotProduct(v0x, v0y, v0z, v0x, v0y, v0z);
  		// dot01 = dot(v0, v1)
  		double dot01 = FunctionalVector3D.dotProduct(v0x, v0y, v0z, v1x, v1y, v1z);
  		// dot02 = dot(v0, v2)
  		double dot02 = FunctionalVector3D.dotProduct(v0x, v0y, v0z, v2x, v2y, v2z);
  		// dot11 = dot(v1, v1)
  		double dot11 = FunctionalVector3D.dotProduct(v1x, v1y, v1z, v1x, v1y, v1z);
  		// dot12 = dot(v1, v2)
  		double dot12 = FunctionalVector3D.dotProduct(v1x, v1y, v1z, v2x, v2y, v2z);
  
  		//
  		// Compute barycentric coordinates
  		//
  		double invDenom = 1. / (dot00 * dot11 - dot01 * dot01);
  		double u = (dot11 * dot02 - dot01 * dot12) * invDenom;
  		double v = (dot00 * dot12 - dot01 * dot02) * invDenom;
  
  		// Check if point is in triangle
  		if ((MathUtil.compareEpsilon(u, 0., epsilon) >= 0)
  				&& (MathUtil.compareEpsilon(v, 0., epsilon) >= 0) 
  				&& (MathUtil.compareEpsilon(u + v, 1., epsilon) <= 0)) {
  			if (forceCoplanar) {
  				// Triangle's plane equation:
  				// nx = ay * (bz - cz) + by * (cz - az) + cy * (az - bz)
  				// ny = az * (bx - cx) + bz * (cx - ax) + cz * (ax - bx)
  				// nz = ax * (by - cy) + bx * (cy - ay) + cx * (ay - by)
  				// d = - (nx * ax + ny * ay + nz * az)
  
  				// Reuse the dot* variables to prevent memory allocation
  				dot00 = ay * (bz - cz) + by * v0z - cy * v1z;
  				dot01 = az * (bx - cx) + bz * v0x - cz * v1x;
  				dot02 = ax * (by - cy) + bx * v0y - cx * v1y;
  				dot11 = - (dot00 * ax + dot01 * ay + dot02 * az);
  				dot12 = dot00 * px + dot01 * py + dot02 * pz + dot11;
  
  				return MathUtil.isEpsilonZero(dot12, epsilon);
  			}
  			return true;
  		}
  		return false;
  	}