Shader实验室在无需三角函数的极简逆向运动学中介绍过如何避免使用三角函数来计算相互约束的逆向运动，它的实现方式比较特殊，所以本文介绍在求解方式上采用比较普遍的方式，即通过三角函数进行求解。
一般来讲，在约束的方式上会采用以下两种方案：
(绿色部分为前臂)
- 前臂向左弯折。
- 前臂向右弯折。
前臂向左弯折
下图为该情况的示意图：
已知a1与a2两骨骼的长度与目标点E，及原点O(这里的原点不一定为(0, 0)点，只是为了便于说明) 到目标点E的距离r，我们要求取的是角度q1与q2。
过OE作垂直于X轴的直角三角形，先求取角度O1：
通过如下的方式可求取O1的角度值：
$$
O_1=tan^{-1}(\frac{Ey-Oy}{Ex-Ox})
$$
接着，我们求取O2，与α的角度值，我们通过连接所有关节点，这样形成了另一个小的三角形：
根据三角形的余弦定理，可求取α与角O2的值：
因为：
$$
r²=a_1²+a_2²-2a_1a_2cos(α)
$$
所以：
$$
α=cos^{-1}(\frac{a_1^2+a_2^2-r^2}{2a_1a_2})
$$
同理可证：
$$
O_2=cos^{-1}(\frac{r^2+a_1^2-a_2^2}{2ra_1})
$$
所以此时可得到∠q2的值为：
$$
q_2=PI-α
$$
此时∠q1的值为：
$$
q_1=O_1-O_2
$$
前臂向右弯折
下图为该情况的示意图：
已知条件与第一种情形一致，求取q1与q2的值。
但是现在与第一种情形不同的地方在于q1的角度是∠O1与∠O2的和：
即：
$$
q_1=O_1+O_2
$$
q2求取的方式也与第一种方式一致，需要注意的是q2角度的值是负值。
因为q2的取值是负值，如果中间关节用C表示，所以也可以直接计算∠OCE的值即可。不需要再通过PI-∠OCE来求取q2的值。但是在游戏引擎实现过程，需要把CE与CO重合时看成0°
最后效果：
ShaderToy源码：
#define STEP 200.
#define STEP 200.
#define R 0.1
float sdCylinder(vec4 p, float r, float h) {
return max(length(p.xz)-r, abs(p.y)-h);
}
float sdSphere(vec4 p, float r) {
return length(p.xyz)-r;
}
float sdBox(vec4 p, vec3 b) {
vec3 q = abs(p.xyz) - b;
return length(max(q,0.0)) + min(max(q.x,max(q.y,q.z)),0.0);
}
// 创建平移矩阵（沿XYZ轴移动）
mat4 translate(float tx, float ty, float tz) {
return mat4(
1, 0, 0, 0, // 第一列
0, 1, 0, 0, // 第二列
0, 0, 1, 0, // 第三列
tx, ty, tz, 1 // 第四列（平移分量）
);
}
// 向量参数版本
mat4 translate(vec3 t) {
return translate(t.x, t.y, t.z);
}
// 创建均匀缩放矩阵
mat4 scale(float s) {
return mat4(
s, 0, 0, 0, // 第一列
0, s, 0, 0, // 第二列
0, 0, s, 0, // 第三列
0, 0, 0, 1 // 第四列
);
}
// 非均匀缩放（各轴独立缩放）
mat4 scale(vec3 s) {
return mat4(
s.x, 0, 0, 0,
0, s.y, 0, 0,
0, 0, s.z, 0,
0, 0, 0, 1
);
}
mat4 rotateX(float theta) {
float c = cos(theta);
float s = sin(theta);
return mat4(
1, 0, 0, 0,
0, c, s, 0, // Y和Z轴的变化
0,-s, c, 0,
0, 0, 0, 1
);
}
mat4 rotateY(float theta) {
float c = cos(theta);
float s = sin(theta);
return mat4(
c, 0,-s, 0, // X和Z轴的变化
0, 1, 0, 0,
s, 0, c, 0,
0, 0, 0, 1
);
}
mat4 rotateZ(float theta) {
float c = cos(theta);
float s = sin(theta);
return mat4(
c, s, 0, 0, // X和Y轴的变化
-s, c, 0, 0,
0, 0, 1, 0,
0, 0, 0, 1
);
}
float smin( float a, float b, float k ) {
float h = clamp( 0.5+0.5*(b-a)/k, 0., 1. );
return mix( b, a, h ) - k*h*(1.0-h);
}
float smax(float a, float b, float k) {
return smin(a, b, -k);
}
float hash(float x) {
return fract(sin(x * 13452.) * 1000.);
}
float hash2(vec2 v) {
return fract(sin(dot(v, vec2(64.24232, 87.873)))*10232.);
}
float hash3(vec3 v) {
return fract(sin(dot(v, vec3(64.24232, 87.873, 76.635)))*10232.);
}
float sdCapsule( vec3 p, vec3 a, vec3 b, float r )
{
vec3 pa = p - a, ba = b - a;
float h = clamp( dot(pa,ba)/dot(ba,ba), 0.0, 1.0 );
return length( pa - ba*h ) - r;
}
float valueNoise(vec2 v) {
vec2 f = floor(v);
vec2 r = smoothstep(vec2(0.), vec2(1.), fract(v));
float lb = hash2(f);
float rb = hash2(f+vec2(1.,0.));
float lt = hash2(f+vec2(0.,1.));
float rt = hash2(f+vec2(1.,1.));
float b_line = mix(lb, rb, r.x);
float t_line = mix(lt, rt, r.x);
return mix(b_line, t_line, r.y);
}
float valueNoise3(vec3 v) {
vec3 f = floor(v);
vec3 r = smoothstep(vec3(0.), vec3(1.), fract(v));
return mix(
mix(mix(hash3(f), hash3(f+vec3(1.,0.,0.)), r.x),
mix(hash3(f + vec3(0.,1.,0.)), hash3(f + vec3(1.,1.,0.)), r.x),
r.y),
mix(mix(hash3(f+vec3(0.,0.,1.)), hash3(f+vec3(1.,0.,1.)), r.x),
mix(hash3(f + vec3(0.,1.,1.)), hash3(f + vec3(1.,1.,1.)), r.x),
r.y),
r.z
);
}
void calcIK(inout vec4 A, inout vec4 B, inout vec4 C,
inout vec4 jaw1, inout vec4 jaw2,vec4 target) {
vec4 lA = A, lB = B, lC = C;
mat4 combineMat;
mat4 wAM = translate(vec3(A.x, A.y, A.z));
vec4 wB = wAM * B;
mat4 wBM = wAM * translate(vec3(B.xyz));
vec4 wC = wBM * C;
float LA = length(wB.xyz - A.xyz);
float LB = length(wC.xyz - wB.xyz);
vec3 AT = target.xyz - A.xyz;
float LT = length(AT);
float diffTy = target.y - A.y;
float diffTx = target.x - A.x;
float diffTz = target.z - A.z;
float zRotAngle = atan(diffTy, sqrt(diffTx*diffTx + diffTz*diffTz));
float yRotAngle = atan(diffTz, diffTx);
if((LA+LB) < LT) {
wAM *= rotateY(-yRotAngle) * rotateZ(zRotAngle);
combineMat = wAM * translate(lB.xyz);
C = combineMat * lC;
float rotateVal = iTime*4.;
mat4 jawMat = combineMat * translate(lC.xyz) * rotateX(rotateVal);
jaw1 = jawMat * jaw1;
jaw2 = jawMat * jaw2;
B = wAM * lB;
} else {
float LA2 = LA*LA;
float LT2 = LT*LT;
float LB2 = LB*LB;
float z2RotAngle = acos((LA2 + LT2-LB2)/(2.*LA*LT));
float zBRotAngle = -(PI - acos((LA2 + LB2-LT2)/(2.*LA*LB)));
wAM *= rotateY(-yRotAngle) * rotateZ(zRotAngle+z2RotAngle);
combineMat = wAM * translate(lB.xyz) * rotateZ(zBRotAngle);
C = combineMat * lC;
float rotateVal = iTime*4.;
mat4 jawMat = combineMat * translate(lC.xyz) * rotateX(rotateVal);
jaw1 = jawMat * jaw1;
jaw2 = jawMat * jaw2;
B = wAM * lB;
}
}
vec2 dist(vec3 p) {
float mat = 0.;
vec4 p_4 = vec4(p,1.);
// joint worldPos
vec4 A = vec4(0., -.3, 0., 1.);
// joint localPos
vec4 B = vec4(1., 0., 0., 1.);
vec4 C = vec4(.8, 0., 0., 1.);
// jaw
vec4 jaw1 = vec4(0., 0., .1, 1.);
vec4 jaw2 = vec4(0.25, 0., .1, 1.);
// world pos
float yPos =mix(hash(floor(iTime)), hash(floor(iTime))+1., fract(iTime));
float h = 0.1;
h =mix(h, 0.8, smoothstep(2., 3., mod(iTime, 6.)) - smoothstep(3., 4., mod(iTime, 6.)));
vec4 tar = vec4(sin(iTime)*1.5, abs(sin(iTime*5.))*h, (cos(iTime))*0.6, 1.);
float d = 0.;
float target = sdSphere(vec4(p_4.xyz - tar.xyz, 1.), 0.1);
calcIK(A, B, C, jaw1, jaw2, tar);
float armT = sdCapsule(p, A.xyz, B.xyz, R);
float armB = sdCapsule(p, B.xyz, C.xyz, R);
d = min(target, min(armT, armB));
vec2 arm = vec2(d, 2.);
float t = iTime;
// 底座
d = sdCylinder(vec4(p - vec3(0,-0.7,0), 1.), 1.3, 0.01);
// 立柱
d = smin(d, sdCylinder(vec4(p - vec3(0,-0.8,0), 1.), 0.2, 0.4), 0.1);
// 第一关节
vec4 p1 = (vec4(p - A.xyz, 1.));
d = smin(d, sdSphere(p1, 0.15), 0.1);
// 大臂
d = smin(d, sdCapsule(p.xyz, A.xyz, B.xyz, 0.1), 0.1);
// 第二关节
vec4 p3 = vec4(p.xyz - B.xyz, 1.);
d = smin(d, sdSphere(p3, 0.12), 0.03);
// 小臂
d = smin(d, sdCapsule(p.xyz, B.xyz, C.xyz, 0.1), 0.03);
// 夹爪
vec3 dirCB = normalize(B.xyz-C.xyz);
d = smin(d, sdCapsule(p.xyz, C.xyz+dirCB*0.07, C.xyz, 0.14), 0.01);
d = smin(d, sdCapsule(p.xyz, jaw1.xyz, jaw2.xyz, 0.04), 0.03);
vec3 offset = normalize(C.xyz - jaw1.xyz)*0.1;
vec3 jaw1_ = C.xyz + offset;
vec3 jaw2_ = jaw2.xyz + 2. * offset;
d = smin(d, sdCapsule(p.xyz, jaw1_, jaw2_, 0.04), 0.03);
float targetSphere = sdSphere(vec4(p.xyz-jaw2_ + offset, 1.), 0.06);
if(targetSphere < d) {
mat = 1.;
}
d = min(d, targetSphere);
return vec2(d, mat);
}
vec3 calNormal(vec3 p) {
vec2 k = vec2(0.001, 0.);
return normalize(vec3(
dist(p + k.xyy).x - dist(p - k.xyy).x,
dist(p + k.yxy).x - dist(p - k.yxy).x,
dist(p + k.yyx).x - dist(p - k.yyx).x
));
}
vec2 raymarch(vec3 ro, vec3 rd) {
float t = 0.;
float y = -1.;
for(float i = 0.; i<STEP; i++) {
vec2 d = dist(ro + rd*t);
if(d.x < .001) {
y = d.y;
break;
}
if(t > STEP) {
t = -1.;
break;
}
t += d.x;
}
return vec2(t, y);
}
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
// Normalized pixel coordinates (from -1 to 1)
vec2 uv = 2. * (fragCoord - 0.5*iResolution.xy) / iResolution.y;
vec3 skyCol = vec3(0.302,0.298,0.298);
vec3 col = skyCol;
vec3 ro = vec3(0.,-.2,1.9);
vec3 rd = normalize(vec3(uv.xy, -1.));
vec2 t = raymarch(ro, rd);
if(t.x > 0.) {
vec3 mate = vec3(0.2);
vec3 sunCol = vec3(1.);
vec3 p = ro + rd*t.x;
vec3 n = calNormal(p);
if(t.y < .5) {
vec3 sunDir = normalize(vec3(0.8,0.4,0.2));
float sunDiff = max(0., dot(n, sunDir));
float skyDiff = clamp(0.5 + 0.5*dot(vec3(0.,1.,0.), n), 0., 1.);
float sunShadow = step(raymarch(p + n*0.01, sunDir).x, 0.);
col = mate * vec3(8.0,5.0,3.0) * sunDiff * sunShadow;
col += mate * skyDiff * vec3(0.2, 0.6, 1.);
} else if(t.y < 1.5) {
vec3 viewDir = normalize(ro - p);
col = pow(1.-dot(viewDir, n)+.8, 5.) * vec3(0.2, 0.6, 1.);
}
}
col = pow(col, vec3(0.4545));
// Output to screen
fragColor = vec4(col,1.0);
}
如果想了解3D逆向运动学或多关节的更多知识可查看：
