本文效果如下:
// The MIT License
// Copyright © 2013 Inigo Quilez
// Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
// An analytical IK solver for 2 segments, WITHOUT trigonometry - gven two
// base and tip white points at which the bones start and end, the red join
// (yellow) is computed automatically by the solver. You can read the
// derivation here: https://iquilezles.org/articles/simpleik/
// A 3D version of this is used in the Insect shader https://www.shadertoy.com/view/Mss3zM
//
// Move the mouse to see the behavior of the solver.
#if 1
// https://iquilezles.org/articles/simpleik/
vec2 solve( vec2 p, float r1, float r2 )
{
vec2 q = p*(0.5+0.5*(r1*r1-r2*r2)/dot(p,p));
float s = r1*r1/dot(q,q)-1.0;
if( s<0.0 ) return vec2(-100.0); // no solution
return q + vec2(-q.y,q.x)*sqrt(s);
}
#else
// alternative, single division
vec2 solve( vec2 p, float r1, float r2 )
{
float h2 = dot(p,p);
float w = h2 + r1*r1 - r2*r2;
float s = 4.0*r1*r1*h2 - w*w;
if( s<0.0 ) return vec2(-100.0); // no solution
return (w*p + vec2(-p.y,p.x)*sqrt(s)) * 0.5/h2;
}
#endif
// https://iquilezles.org/articles/distfunctions2d/
float sdSegment( vec2 a, vec2 b, vec2 p )
{
vec2 pa = p - a;
vec2 ba = b - a;
float h = clamp( dot(pa,ba)/dot(ba,ba), 0.0, 1.0 );
return length( pa - ba*h );
}
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
// pixel setup
vec2 uv = (2.0*fragCoord-iResolution.xy) / iResolution.y;
vec2 mo = (2.0*iMouse.xy-iResolution.xy) / iResolution.y;
if( iMouse.z<=0.0001 ) { mo = vec2(0.8+0.3*sin(iTime),0.2+0.2*sin(iTime*1.31)+0.2*sin(iTime*2.11)); }
// compute IK
const float l1 = 0.8; // length of first segment
const float l2 = 0.4; // length of second segment
vec2 a = vec2(0.0,0.0); // location of base
vec2 c = mo; // location of tip
vec2 b = solve( c, l1, l2 ); // compute location of join
// draw background
vec3 col = vec3(0.3)*(1.0-0.2*length(uv));
// plot circles
col = mix( col, vec3(0.5,0.5,0.5), 1.0-smoothstep( 0.0, 0.007, abs(length( a - uv )-l1) ) );
col = mix( col, vec3(0.5,0.5,0.5), 1.0-smoothstep( 0.0, 0.007, abs(length( c - uv )-l2) ) );
// plot segments and join
if( b.x>-99.0 )
{
float f = min( sdSegment( a, b, uv ), sdSegment( b, c, uv ) ) - 0.002;
col = mix( col, vec3(1.0,1.0,1.0), 1.0-smoothstep( 0.000, 0.005, f ) );
col = mix( col, vec3(1.0,1.0,0.0), 1.0-smoothstep( 0.035, 0.040, length( b - uv ) ) );
}
// plot base and tip
col = mix( col, vec3(1.0,1.0,1.0), 1.0-smoothstep( 0.035, 0.040, length( a - uv ) ) );
col = mix( col, vec3(1.0,1.0,1.0), 1.0-smoothstep( 0.035, 0.040, length( c - uv ) ) );
// dither to remove banding in the background
col += fract(sin(fragCoord.x*vec3(13,17,11)+fragCoord.y*vec3(1,7,5))*158.391832)/255.0;
fragColor = vec4( col, 1.0 );
}
IK代表"逆向运动学
"(Inverse Kinematics)
,这是一个复杂的课题。但本文讨论的是一种特殊的简化模型,实际上是可能实现的最极端简化形式,这种简化使得该问题能够获得解析解。在这种简化下,问题可总结为:空间中存在一个固定点,两个已知长度(r1和r2)的骨骼段通过连接点x串联,其中第一段长度r1连接固定点与连接点x,第二段长度r2连接x与目标位置p。我们的目标是仅根据已知的目标点p坐标以及长度r1和r2,确定连接点x在空间中的位置(即其坐标)。
遗憾的是,若您查阅相关解决方案资料,会发现人们普遍使用三角函数来解决这一问题(即确定x的位置)。我认为这种现状存在多重弊端,其中一个关键原因在于:我们完全可以通过向量运算与几何原理,以更加简洁优雅的方式实现求解。尽管我此前未曾见过类似推导,但若他人也以相同方式得出此解,我并不会感到意外。尽管如此,我仍决定撰写本文专门阐述这种逆向运动学求解方法——因为若不主动探寻,您很可能永远无法偶然发现它。方法如下:
解决方案
我们首要任务是精确定位求解点x在空间中的位置。为简化示意图,我将移除原图中的坐标轴,转而绘制两个关键圆:根据几何约束,求解点x必须满足以下条件——
- 距离原点r₁个单位,故必定位于以原点为圆心、r₁为半径的圆周上
- 同时距离目标点p为 r₂个单位,故亦需处于以p为中心、r₂为半径的圆周上
这两个圆的交点即为潜在求解点x的可能位置。
通过这种几何视角,可以直观发现:由于两圆相交于两个不同位置,解点x实际上存在两个潜在解。为精确求解这些交点,我们需将几何图示转化为数学方程——即分别为两个圆建立方程:
-
原点约束圆
$$
|x|^2=r1^2
$$(描述x距离原点为r₁的所有可能位置)
-
目标点约束圆
$$
|x-p|^2=r2^2
$$(描述x距离目标点p为r₂的所有可能位置)
现在我们可以展开第二个方程中的二项式得到:
$$
|x|^2-2(x·p)+|p|^2=r2^2
$$
然后我们可以把|x|²从第一个方程(原点约束圆)替换到展开后的方程中:
$$
r1² - 2(x⋅p) + |p|² = r2²
$$
然后重新排列后就得到了:
$$
(x⋅p) = ( r1² - r2² + |p|²)/2
$$
这是我们的第一组约束方程。此时方程组仍包含两个变量(即x点的二维坐标),因此需从几何关系中提取更多信息方能求解。关键几何特性在于:
两解点连线必定与原点-p点连线相互垂直。换言之,存在一个位于原点-p连线上的特殊点q,该点满足:
- 是两解点的精确中点
- 可通过将x点投影到原点到p点连线上求得
那么,我们通过点积来计算在q点的投影:
$$
q = p (x⋅p) / |p|²
$$
现在我们可以把前面计算的表达式(x·p)代入上式得到公式1:
$$
q = p ½ + ½(r1² - r2²) / |p|²
$$
有点接近答案了。看一下这个图,我们可以把x表示为下式,得到公式2:
$$
x = q ± s q⊥(公式2)
$$
对于某个标量s(其中q⊥表示q的垂直分量),s本质上表征了从点q沿q垂直方向到x的垂直偏移量,其度量单位以|q|为基准。值得庆幸的是,通过将x的这一几何描述代入最初的约束方程(即半径为r₁的圆方程),我们即可轻松求得s的值:
$$
|x|² = r1²
$$
得到了:
$$
|q ± s q⊥|² = r1²
$$
然后二项式展开:
$$
|q|² ± 2s(q⋅q⊥) + s²|q⊥|² = r1²
$$
并继续注意到,q和q⊥之间的点积为零,因为它们是垂直的。我们还注意到,根据定义,q⊥的长度与q的长度相同
(没理解,指的特定条件
?),因此我们可以将上述方程简化为:
$$
|q|² (1 + s²) = r1²
$$
最后解出s,得到公式3:
$$
s = ±\sqrt{(r1²/|q|² - 1)}(公式3)
$$
实际上要计算s很简单,因为我们已知原点约束圆的公式以及q点,根据勾股定理:
$$
s^2=|x|^2-|q|^2
$$
所以s就等于:
$$
s=±\sqrt{(r1²-|q|^2)}
$$
由于|q|²等于dot(q, q),所以 只要结合公式1
就可以求得s:
$$
s=±\sqrt{(r1²-dot(q, q))}(公式2)
$$
这里的平方根的两个解代表了点x的两个可能位置。当然,如果平方根的参数是负数,则意味着这两个圆不相交,我们的问题没有解,即我们移动点p到一个程度,使得使用我们现有的长度为r1和r2的两段无法将其连接到原点。在这种情况下,我们可以向调用者发出信号表示不存在解,或者我们可以强制让s为非负值,这相当于拉伸r1和r2的长度。当拉伸是最小的,并且仅由于我们浮点计算中的微小精度问题而发生时,这应该是一个可行的解决方案。
代码
代码是对上述公式1、2和3的简单实现:
vec2 solve( vec2 p, float r1, float r2 )
{
vec2 q = p*(0.5+0.5*(r1*r1-r2*r2)/dot(p,p));
float s = r1*r1/dot(q,q)-1.0;
s = max(s,0.0);
//if( s<0.0 ) return vec2(-inf); // or we could instead signal no solution
return q + vec2(-q.y,q.x)*sqrt(s);
}
在上述实现中,我选择了正平方根,但你可能需要灵活处理,并选择对你的特定配置最有意义的那个。
无论如何,我们可以通过一些代数方法来稍微优化这段代码,使得我们可以用单次除法重写上面的代码。
vec2 solve( vec2 p, float r1, float r2 )
{
float h = dot(p,p);
float w = h + r1*r1 - r2*r2;
float s = max(4.0*r1*r1*h - w*w,0.0);
return (w*p + vec2(-p.y,p.x)*sqrt(s)) * 0.5/h;
}
在三维空间计算中我们新增IK骨骼运动的方向,结合上述改进后的公式1与公式2,得到
:
vec3 solve( vec3 p, float r1, float r2, vec3 dir )
{
vec3 q = p*( 0.5 + 0.5*(r1*r1-r2*r2)/dot(p,p) );
float s = r1*r1 - dot(q,q);
s = max( s, 0.0 );
q += sqrt(s)*normalize(cross(p,dir));
return q;
}
最后,这里是在Shadertoy中运行的代码供参考,带有一些鼠标交互功能,所以你可以看到系统如何响应不同的参数配置:
// The MIT License
// Copyright © 2013 Inigo Quilez
// Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
// An analytical IK solver for 2 segments, WITHOUT trigonometry - gven two
// base and tip white points at which the bones start and end, the red join
// (yellow) is computed automatically by the solver. You can read the
// derivation here: https://iquilezles.org/articles/simpleik/
// A 3D version of this is used in the Insect shader https://www.shadertoy.com/view/Mss3zM
//
// Move the mouse to see the behavior of the solver.
#if 1
// https://iquilezles.org/articles/simpleik/
vec2 solve( vec2 p, float r1, float r2 )
{
vec2 q = p*(0.5+0.5*(r1*r1-r2*r2)/dot(p,p));
float s = r1*r1/dot(q,q)-1.0;
if( s<0.0 ) return vec2(-100.0); // no solution
return q + vec2(-q.y,q.x)*sqrt(s);
}
#else
// alternative, single division
vec2 solve( vec2 p, float r1, float r2 )
{
float h2 = dot(p,p);
float w = h2 + r1*r1 - r2*r2;
float s = 4.0*r1*r1*h2 - w*w;
if( s<0.0 ) return vec2(-100.0); // no solution
return (w*p + vec2(-p.y,p.x)*sqrt(s)) * 0.5/h2;
}
#endif
// https://iquilezles.org/articles/distfunctions2d/
float sdSegment( vec2 a, vec2 b, vec2 p )
{
vec2 pa = p - a;
vec2 ba = b - a;
float h = clamp( dot(pa,ba)/dot(ba,ba), 0.0, 1.0 );
return length( pa - ba*h );
}
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
// pixel setup
vec2 uv = (2.0*fragCoord-iResolution.xy) / iResolution.y;
vec2 mo = (2.0*iMouse.xy-iResolution.xy) / iResolution.y;
if( iMouse.z<=0.0001 ) { mo = vec2(0.8+0.3*sin(iTime),0.2+0.2*sin(iTime*1.31)+0.2*sin(iTime*2.11)); }
// compute IK
const float l1 = 0.8; // length of first segment
const float l2 = 0.4; // length of second segment
vec2 a = vec2(0.0,0.0); // location of base
vec2 c = mo; // location of tip
vec2 b = solve( c, l1, l2 ); // compute location of join
// draw background
vec3 col = vec3(0.3)*(1.0-0.2*length(uv));
// plot circles
col = mix( col, vec3(0.5,0.5,0.5), 1.0-smoothstep( 0.0, 0.007, abs(length( a - uv )-l1) ) );
col = mix( col, vec3(0.5,0.5,0.5), 1.0-smoothstep( 0.0, 0.007, abs(length( c - uv )-l2) ) );
// plot segments and join
if( b.x>-99.0 )
{
float f = min( sdSegment( a, b, uv ), sdSegment( b, c, uv ) ) - 0.002;
col = mix( col, vec3(1.0,1.0,1.0), 1.0-smoothstep( 0.000, 0.005, f ) );
col = mix( col, vec3(1.0,1.0,0.0), 1.0-smoothstep( 0.035, 0.040, length( b - uv ) ) );
}
// plot base and tip
col = mix( col, vec3(1.0,1.0,1.0), 1.0-smoothstep( 0.035, 0.040, length( a - uv ) ) );
col = mix( col, vec3(1.0,1.0,1.0), 1.0-smoothstep( 0.035, 0.040, length( c - uv ) ) );
// dither to remove banding in the background
col += fract(sin(fragCoord.x*vec3(13,17,11)+fragCoord.y*vec3(1,7,5))*158.391832)/255.0;
fragColor = vec4( col, 1.0 );
}
结论
如你所见,我们通过几何推理而非三角函数解决了这个问题,而且我们的实现也没有使用任何三角函数。通常情况下,正如我在其他时候表达的那样,如果你发现自己在解决某些二维或三维几何问题时使用了三角函数,那么你可能采取了一种不够优雅和高效的方法。当然,这只是一个观点问题,但现在你已经听到了我的看法😄
无论如何,这种逆向运动学(IK)公式在为关节对象(例如这个ShaderToy例子中的动画生物的腿)做简单动作时非常有用。你可以点击播放来查看它的动画效果,或者点击标题查看源代码:示例链接
原文作者:Inigo Quilez