Simple illustration of the Jacobian determinant as a geometric scaling factor
- Consider a simple quadratic mapping from x to X
- We start with an interval in x from a to b
- Then find the equivalent interval in X from the mapped points
- Note that the ratio of the lengths is the gradient d(x^2)/dx = 2x
- In two dimensions, we can consider transforming a unit square
- Consider a linear transformation of the coordinates (with a matrix)
- The parallelogram can be divided into squares and triangles
- We can also check this numerically using a MATLAB function
- Importantly, the area is the determinant of the matrix T
- The extension to three dimensions is analogous
- For a nonlinear transformation think back to the slope of the quadratic
- The matrix elements are the transformation's partial derivatives
Mathematically, this is an example of the change of variable theorem for integrals of multidimensional functions .
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.
(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')
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')
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')
The gradient is evaluated in the middle of our chosen interval.
grad = 2*x; grad5 = grad(x == m) 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.
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')
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')
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
A = p*s - q*r A = polyarea(COORDS(1, :), COORDS(2, :))
A = 0.9300 A = 0.9300
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 unit cube maps to a parallelepiped, whose volume is given by the determinant of the matrix that represents the linear transformation .
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.
[dX1/dx1 dX1/dx2 dX1/dx3] J = [dX2/dx1 dX2/dx2 dX2/dx3] [dX3/dx1 dX3/dx2 dX3/dx3]
-  http://en.wikipedia.org/wiki/Integration_by_substitution
-  http://en.wikibooks.org/wiki/SPM/VBM#.22Modulation.22
-  http://www.fil.ion.ucl.ac.uk/spm/doc/books/hbf2/pdfs/Ch6.pdf
-  http://www.math.uah.edu/stat/dist/Transformations.xhtml
-  http://www.fil.ion.ucl.ac.uk/~wpenny/publications/pdfs.ps
-  http://en.wikipedia.org/wiki/Parallelepiped
Copyright 2010 Ged Ridgway