Fix: Pole Vector chooses correct bend axis

Now, we find out where the "knee" is in the IK chain and check which of
the base-bone's axes it is in-line with. Then we calculate the poll
vector based around the perpendicular axis.
My simple tests worked perfectly. Need to test this more since the
math is confusing.
Also added some math utilities to utilities.py. Bisect isn't used
but since I went to the trouble to find it and it is useful, I'm
keeping it in, in case I need it later.

link_containers.py

@@ -1548,7 +1548,7 @@ class LinkInverseKinematics:
         print(wrapGreen("Creating ")+wrapOrange("Inverse Kinematics")+
              wrapGreen(" Constraint for bone: ") +
              wrapOrange(self.GetxForm().bGetObject().name))
-        myOb = self.GetxForm().bGetObject()
+        ik_bone = self.GetxForm().bGetObject()
         c = self.GetxForm().bGetObject().constraints.new('IK')
         get_target_and_subtarget(self, c)
         get_target_and_subtarget(self, c, input_name = 'Pole Target')
@@ -1559,22 +1559,40 @@ class LinkInverseKinematics:
         c.chain_count = 1 # so that, if there are errors, this doesn't print a whole bunch of circular dependency crap from having infinite chain length
         if (c.pole_target): # Calculate the pole angle, the user shouldn't have to.
             pole_object = c.pole_target
-            pole_location = pole_object.matrix_world.decompose()[0]
+            assert pole_object == c.target, f"Error with {self}: Pole Target must be bone within the same Armature as IK Bone  -- for now."
+            pole_location = None
             if (c.pole_subtarget):
                 pole_object = c.pole_target.pose.bones[c.pole_subtarget]
-                pole_location += pole_object.bone.head_local
+                pole_location = pole_object.bone.head_local
+            else: #TODO this is a dumb limitation but I don't want to convert to the armature's local space so that some idiot can rig in a stupid way
+                raise RuntimeError(f"Error with {self}: Pole Target must be bones within the same Armature as IK Bone -- for now.")
             #HACK HACK
-            handle_location = myOb.bone.tail_local if (self.evaluate_input("Use Tail")) else myOb.bone.head_local
+            handle_location = ik_bone.bone.tail_local if (self.evaluate_input("Use Tail")) else ik_bone.bone.head_local
             counter = 0
-            parent = myOb
-            base_bone = myOb
+            parent = ik_bone
+            base_bone = ik_bone
             while (parent is not None):
+                counter+=1
                 if ((self.evaluate_input("Chain Length") != 0) and (counter > self.evaluate_input("Chain Length"))):
                     break
                 base_bone = parent
                 parent = parent.parent
-                counter+=1
-            head_location = base_bone.bone.head_local
+
+            def get_main_axis(bone, knee_location):
+                # To decide whether the IK mainly bends around the x or z axis....
+                x_axis = bone.matrix_local.to_3x3() @ Vector((1,0,0))
+                y_axis = bone.matrix_local.to_3x3() @ Vector((0,1,0))
+                z_axis = bone.matrix_local.to_3x3() @ Vector((0,0,1))
+                # project the knee location onto the plane of the bone.
+                from .utilities import project_point_to_plane
+                planar_projection = project_point_to_plane(knee_location, bone.head_local, y_axis)
+                # and get the dot between the X and Z axes to find which one the knee is displaced on.
+                x_dot = x_axis.dot(planar_projection) # whichever axis' dot-product is closer to zero 
+                z_dot = z_axis.dot(planar_projection) #  with the base_bone's axis is in-line with it.
+                prWhite(bone.name, z_dot, x_dot)
+                # knee is in-line with this axis vector, the bend is happening on the perpendicular axis.
+                if abs(z_dot) < abs(x_dot): return x_axis # so we return X if Z is in-line with the knee
+                else: return z_axis                       # and visa versa
 
             # modified from https://blender.stackexchange.com/questions/19754/how-to-set-calculate-pole-angle-of-ik-constraint-so-the-chain-does-not-move
             from mathutils import Vector
@@ -1586,16 +1604,26 @@ class LinkInverseKinematics:
                     angle = -angle
                 return angle
 
-            def get_pole_angle(base_bone, ik_bone, pole_location):
+            def get_pole_angle(base_bone, ik_bone, pole_location, main_axis):
                 pole_normal = (ik_bone.bone.tail_local - base_bone.bone.head_local).cross(pole_location - base_bone.bone.head_local)
                 projected_pole_axis = pole_normal.cross(base_bone.bone.tail_local - base_bone.bone.head_local)
-                x_axis= base_bone.bone.matrix_local.to_3x3() @ Vector((1,0,0))
-                return signed_angle(x_axis, projected_pole_axis, base_bone.bone.tail_local - base_bone.bone.head_local)
-
-            pole_angle_in_radians = get_pole_angle(base_bone,
-                                                myOb,
-                                                pole_location)
+                # note that this normal-axis is the y-axis but flipped
+                return signed_angle(main_axis, projected_pole_axis, base_bone.bone.tail_local - base_bone.bone.head_local)
+
+            if self.evaluate_input("Use Tail") == True:
+                main_axis = get_main_axis(ik_bone.bone, ik_bone.bone.tail_local)
+                # pole angle to the PV:
+                pole_angle_in_radians = get_pole_angle(base_bone, ik_bone, pole_location, main_axis)
+            elif ik_bone.bone.parent:
+                main_axis = get_main_axis(ik_bone.bone.parent, ik_bone.bone.tail_local)
+                pole_angle_in_radians = get_pole_angle(base_bone, ik_bone, pole_location, main_axis)
+            else: # the bone is not using "Use Tail" and it has no parent -- meaningless.
+                pole_angle_in_radians = 0
+            
             c.pole_angle = pole_angle_in_radians
+
+            # TODO: the pole target should be a bone in a well-designed rig, but I don't want to force this, so....
+            #   in future, calculate all this in world-space so we can use other objects as the pole.
         
         props_sockets = {
         'chain_count'   : ("Chain Length", 1),

utilities.py

@@ -690,7 +690,8 @@ def SugiyamaGraph(tree, iterations):
                     n.location.x -= next_node.width*1.5
         
 
-
+def project_point_to_plane(point, origin, normal):
+    return point - normal.dot(point- origin)*normal
 
 ##################################################################################################
 # stuff I should probably refactor!!
@@ -1185,3 +1186,66 @@ def data_from_ribbon_mesh(m, factorsList, mat, ribbons = None, fReport = None):
             normals.append( outNorm )
         ret.append( (points, widths, normals) )
     return ret # this is a list of tuples containing three lists
+
+
+
+
+#This bisection search is generic, and it searches based on the
+# magnitude of the error, rather than the sign.
+# If the sign of the error is meaningful, a simpler function
+# can be used.
+def do_bisect_search_by_magnitude(
+        owner, 
+        attribute,
+        index = None,
+        test_function = None,
+        modify = None,
+        max_iterations = 10000,
+        threshold = 0.0001,
+        thresh2   = 0.0005,
+        context = None,
+        update_dg = None,
+    ):
+    from math import floor
+    i = 0; best_so_far = 0; best = float('inf')
+    min = 0; center = max_iterations//2; max = max_iterations
+    # enforce getting the absolute value, in case the function has sign information
+    # The sign may be useful in a sign-aware bisect search, but this one is more robust!
+    test = lambda : abs(test_function(owner, attribute, index, context = context,))
+    while (i <= max_iterations):
+        upper = (max - ((max-center))//2)
+        modify(owner, attribute, index, upper, context = context); error1 = test()
+        lower = (center - ((center-min))//2)
+        modify(owner, attribute, index, lower, context = context); error2 = test()
+        if (error1 < error2):
+            min = center
+            center, check = upper, upper
+            error = error1
+        else:
+            max = center
+            center, check = lower, lower
+            error = error2
+        if (error <= threshold) or (min == max-1):
+            break
+        if (error < thresh2):
+            j = min
+            while (j < max):
+                modify(owner, attribute, index, j * 1/max_iterations, context = context)
+                error = test()
+                if (error < best):
+                    best_so_far = j; best = error
+                if (error <= threshold):
+                    break
+                j+=1
+            else: # loop has completed without finding a solution
+                i = best_so_far; error = test()
+                modify(owner, attribute, index, best_so_far, context = context)
+                break
+        if (error < best):
+            best_so_far = check; best = error
+        i+=1
+        if update_dg:
+            update_dg.update()
+    else: # Loop has completed without finding a solution
+        i = best_so_far
+        modify(owner, attribute, best_so_far, context = context); i+=1