The pivot, p, is defined as the point, within the convex hull of the
contacts, which, after rotation and scaling, leaves the transformed
contacts as close to the actual positions as possible.

Let r_i be the starting points and s_i the actual ending points in a
transformation. Let D be the scaling, and R the rotation. Then, minimizing

L(p) = sum_i |D R (r_i - p) + p - s_i|^2 / N

yields the pivot. Let

rm = sum_i r_i / N,
p = rm + u,
q_i = s_i - rm - D R (r_i - rm),

and we get

L(p) = sum_i |(1 - D R) u - q_i|^2 / N.

With

L0 = sum_i norm2(q_i) / N,
T = (1 - D R)' (1 - D R),
m = sum_i q_i / N,

we can write this as

L(p) = L0 + u' T u - 2 m' (1 - D R) u.

To handle the constraint, we can approximate the hull with a circle
centered at rm. If we pick the average radius, P, the constraint becomes

|u| < P.

Relaxing [1] the expression (h >= 0) yields

L(p, h) = L0 + u' T u - 2 m' (1 - D R) u + h (|u|^2 - P^2),

leading to the linear equation

(T + h) u = (1 - D R)' m.

Further,

sm = sum_i s_i / N,
m = sum_i (s_i - rm - D R (r_i - rm)) / N = sm - rm,

thus m is the average displacement. In words, the pivot is the average
position plus a correction depending on the average displacement.

*

Some algebra solves the equation,

D' = D,
[D, R] = 0,
R = S + C,
S' = -S,
C' = C,
R + R' = 2 C,
T = (1 - D R)' (1 - D R) = 1 + D^2 - 2 D C,

which is a simple diagonal scaling operator. With

a = 1 - D C,
b = D S,

we can write this as

T = (1 - DC)^2 + D^2(1 - C^2) = (1 - DC)^2 + D^2 S^2 = a^2 + b^2.

Similarly, we can write

(1 - D R)' = ((a, b), (-b, a)),

and thusly,

u = Q(h) m,

with

Q(h) = ((a, b), (-b, a)) / (a^2 + b^2 + h).

When D R = 1, it follows that a^2 + b^2 = 0, and the relaxation ensures
that u is finite.

*

The drag is found by minimizing

E(d) = sum_i | D R (r_i - p) + p + d - s_i |^2 / N,
E(d) = d^2 + 2 d' ((1 - D R) u - m) + E0,

which leads to the linear equation

d = m - (1 - D R) u.

Explicitly,

d = m - (a ux - b uy, a uy + b ux).

Inserting the expression for u yields, after some algebra,

d = m (1 - (a^2 + b^2) / (a^2 + b^2 + h)).

When h = 0, d = 0, as expected.

When a^2 + b^2 = 0, d = m, also as expected.

For constrained cases, the drag is a fraction of the average displacement.

*

Time to look at measures for the relaxation parameter. Since d depends on
h, we can write the correction u(h) in terms of d instead. After som
algebra,

|u(h)| = (|m| - |d|) / sqrt(a^2 + b^2).

Conversely, d(h) can be written in terms of the constrained u(h) as

d(h) = m (1 - sqrt(a^2 + b^2) |u(h)| / |m|).

Since |u(0)| = |m| / sqrt(a^2 + b^2), we obtain

d(h) = m (1 - |u(h)| / |u(0)|).

*

We can now write down an explicit recipe for determining the pivot (p) and
drag (d), given the transformation parameters a and b.

w = (a mx + b my, a my - b mx).

If |w| = 0, then u = 0. Consequently p = rm, d = m, and we are done. Else,

u = w |m|^2 / |w|^2,

t = P / |u|.

If t >= 1, then p = rm + u, d = 0, and we are done. Else,

p = rm + t u,
d = (1 - t) m.

[1] See Lagrange relaxation