Code covered by the BSD License

# Modulation tutorial

Simple illustration of the Jacobian determinant as a geometric scaling factor

Simple illustration of the Jacobian determinant as a geometric scaling factor

# Simple illustration of the Jacobian determinant as a geometric scaling factor

## Notes

Mathematically, this is an example of the change of variable theorem for integrals of multidimensional functions [1].

It's useful for understanding tensor- and voxel-based morphometry [2,3], which are often of interest to neuroimaging researchers without mathematical backgrounds, which was my motivation for writing this tutorial.

In VBM, use of the Jacobian determinant scaling factor to preserve local probabilistic volumes in spatially normalised (geometrically transformed) tissue segmentation images is known as modulation .

This is also closely related to finding the probability density function of a transformed random variable in terms of the (here assumed monotonic) transformation and the PDF of the original variable [4,5].

We will use an intuitive graphical visualisation.

## Consider a simple quadratic mapping from x to X

(we avoid y to reduce confusion when considering two or three dimensions, also it then seems intuitive to use capitals for other mapped quantities)

x = linspace(0, 3.5, 351); % (chosen to hit x=2.5 exactly)
X = x.^2;
clf; hold on
plot(x, X)
ylim([0 10])
xlabel('x'); ylabel('X')


## We start with an interval in x from a to b

a = 2; b = 3;
l = b - a; % original length
m = (a + b) / 2; % mid-point
plot([a b], [0 0], 'r', 'LineWidth', 3)
plot([a a], [0 a^2], 'k:')
plot([b b], [0 b^2], 'k:')
text(m, 0.5, sprintf('Length l = %g', l), 'HorizontalAlignment', 'Center')


## Then find the equivalent interval in X from the mapped points

A = a^2; B = b^2;
L = B - A;
M = (A + B) / 2; % (Aside: note this is slightly larger than m^2)
plot([a 0], [A A], 'k:')
plot([b 0], [B B], 'k:')
plot([0 0], [A B], 'r', 'LineWidth', 3)
text(0.1, M, sprintf('Length L = %g', L), 'HorizontalAlignment', 'Left')


## Note that the ratio of the lengths is the gradient d(x^2)/dx = 2x

The gradient is evaluated in the middle of our chosen interval.

grad = 2*x;
ratio = L / l

grad5 =

5

ratio =

5



If we moved the (length-2) interval to the right, where the function is steeper, the mapped interval would be longer, and vice versa. If we changed the function so that it was steeper at x = m = 2.5, then the mapped interval would also get longer.

## In two dimensions, we can consider transforming a unit square

coords = [
0 0
1 0
1 1
0 1
]';  % x1 and x2, anti-clockwise around unit square


which we can visualise in the plane (NB the vertical axis is no longer a mapping of the horizontal; we just visualise the original unmapped 2d object first, then we will visualise the mapped one overlaid later).

clf; hold on
f = fill(coords(1,:), coords(2,:), 'b');
axis([-0.5 1.5 -0.5 1.5]); axis equal
xlabel('x1'); ylabel('x2')


## Consider a linear transformation of the coordinates (with a matrix)

p = 1.1; q = 0.3;
r = 0.2; s = 0.9;
T = [p q
r s]
COORDS = T * coords % NB we have no translation, so (0,0) maps to (0,0)

T =

1.1000    0.3000
0.2000    0.9000

COORDS =

0    1.1000    1.4000    0.3000
0    0.2000    1.1000    0.9000



Note that the coordinates map to

[0 p p+q q]
[0 r r+s s]

and observe that the square is rotated (about the origin), stretched and skewed, into a parallelogram, which we now overlay:

F = fill(COORDS(1, :), COORDS(2, :), 'r');
xlabel('x1 or X1'); ylabel('x2 or X2')
text(p, r, '(p, r)', 'HorizontalAlignment', 'Right', ...
'VerticalAlignment', 'Bottom')
text(p+q, r+s, '(p+q, r+s)', 'VerticalAlignment', 'Bottom')
text(q, s, '(q, s)', 'HorizontalAlignment', 'Left', ...
'VerticalAlignment', 'Top')


## The parallelogram can be divided into squares and triangles

This lets us compute the area algebraically in terms of p, q, r and s.

delete(f);
xlabel('X1'); ylabel('X2')
fill([0 p p], [0 0 r], 'g');
fill([p p+q p+q p], [0 0 r r], 'c');
fill([p p+q p+q], [r r r+s], 'm');
fill([q q p+q], [r+s s r+s], 'g');
fill([0 0 q q], [r+s s s r+s], 'c');
fill([0 0 q], [s 0 s], 'm');


The area of the large enclosing square is

The area of each cyan square is q*r; the area of each magenta triangle is q*s/2, and the area of each green triangle is p*r/2, so the red (mapped) area is given by

## We can also check this numerically using a MATLAB function

A = p*s - q*r
A = polyarea(COORDS(1, :), COORDS(2, :))

A =

0.9300

A =

0.9300



## Importantly, the area is the determinant of the matrix T

This also remains true when the shape of the mapped region requires a different construction of squares and triangles to derive the area; it is still p*s - q*r, and this is the determinant of the transformation matrix

A = det(T)

A =

0.9300



## The extension to three dimensions is analogous

The unit cube maps to a parallelepiped, whose volume is given by the determinant of the matrix that represents the linear transformation [6].

## For a nonlinear transformation think back to the slope of the quadratic

If we zoom in enough (mathematically, taking a limit) the quadratic looks like a straight line with its slope equal to the quadratic's gradient at that point. The same is true for higher dimensional nonlinear transformations: if we zoom in enough, they are approximately linear, and their effect on an infinitessimal cube tends to the effect of a matrix multiplication, and the volume scaling factor (between the original and transformed infinitessimal cubes) is the determinant of this matrix.

## The matrix elements are the transformation's partial derivatives

    [dX1/dx1 dX1/dx2 dX1/dx3]
J = [dX2/dx1 dX2/dx2 dX2/dx3]
[dX3/dx1 dX3/dx2 dX3/dx3]