Bundle Adjustment – Part 1: Jacobians

Bundle Adjustment (BA) is a well established technique in Computer Vision typically used as the last step in many feature based 3D reconstruction algorithms to tweak the intrinsic and extrinsic camera parameters as well as the position of the reconstruction points to minimize the reconstruction error between the projected and observed image points¹

In this series of posts, we’ll examine BA in detail. In the first post, we’ll derive expressions for the elements of the Jacobian matrix used in BA, in the second post, we’ll plot some of these Jacobians and consider their properties and in the third and fourth post, we’ll apply BA to simulated data sets and analyze the results.

Assume each camera $j$ is parameterized by a vector $a_j$ and each 3D point $i$ by a vector $b_i$ . Then, Bundle Adjustment involves minimizing the following objective function:

$\min_{a_j,b_i}\sum\limits_{i=1}^n\sum\limits_{j=1}^m d(Q(a_j, b_i), x_{ij})^2$

Here, $Q(a_j, b_i)$ is the predicted projection of point $i$ on the image $j$ and $d(x,y)$ denotes the Euclidean distance between the image points represented by the vectors $x$ and $y$ .

The vector $a_j$ consists of the rotation and translation parameters to transform the world point in the camera coordinate system and the intrinsic parameters of the camera that is the focal length (in the x and y dimensions), coordinates of the center of projection and the skew factor (if any).

Together, these parameters are combined as follows to project a 3D world point on the image:

$\begin{bmatrix}wu\\wv\\w\end{bmatrix} = \begin{bmatrix}f_x & s_k &u_0 \\0 & f_y & v_0 \\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix}R & T \end{bmatrix}\begin{bmatrix}X \\ Y \\ Z \\ 1 \end{bmatrix}$

Henceforth, we’ll assume that our camera is a rectilinear camera with square pixels. Thus, $f_x = f_y$ and $s_k = 1$ .

Bundle Adjustment methods typically employ the Levenberg Marquardt (LM) algorithm to find the minimum of the optimization function. The LM algorithm needs the jacobians i.e., the partial derivatives of the image coordinates wrt the intrinsic and extrinsic parameters of the camera, and the coordinates of the 3D points. During my research on Bundle Adjustment, I found plenty of references that describe the structure of the problem and the techniques used to solve it, but none that show how to compute these jacobians that are needed by any practical implementation. The primary purpose of this post is to show how to compute the jacobians and practical information about some of their properties.

Let’s first consider the number of independent parameters that we need to optimize over. For $N$ points each of which is visible in $M$ images, the number of independent parameters is:

Here, $X, Y, Z$ are the 3D coordinates of the point, $f, u_0, v_0$ are the focal length of the camera and the coordinates of the center of projection respectively (we assume square pixels and zero skew), $w_x, w_y, w_z$ are the parameters of the rotation matrix in the axis-angle representation and $t_x, t_y, t_z$ are the elements of the translation vector between the object and the camera coordinate system. Note that while the rotation matrix has 9 elements, due to the orthogonality and unit norm constraints, only 3 are independent. The axis-angle representation is a convenient format to represent rotations as it only has 3 elements and therefore keeps the number of elements used to represent a rotation to a minimum.

The axis-angle and rotation matrix formats can be inter-converted using the following equations.

$R=e^{\lbrack w\rbrack_x}=I+\frac{sin(\|w\|)}{\|w\|} \lbrack w \rbrack _x + \frac{(1-cos(\|w\|))}{\|w\|^2} (\lbrack w \rbrack _x)^2$

$\|w\| = \arccos(\frac{trace(R)-1}{2}, \frac{w}{\|w\|}=\frac{1}{2sin(\|w\|)}\begin{bmatrix}r_{32}-r_{23}\\r_{13}-r_{31}\\r_{21}-r_{12}\end{bmatrix}$

Here $\lbrack w \rbrack _x$ has the following form:

$\lbrack w \rbrack _x = \begin{bmatrix} 0 & -w_z & w_y \\ w_z & 0 & -w_x \\ -w_y & w_x & 0 \end{bmatrix}$

If we represent the image of a 3D point $P$ by a 2D vector $\bf{x} = (u,v)$ then the jacobian matrix can expressed as two separate matrices $A$ and $B$ that encode the variation of $\bf{x}$ by the camera parameters denoted by the vector $\bf{a}$ and 3D point coordinates denoted by the vector $\bf{b}$ .

$A = \begin{bmatrix}\frac{\partial{u}}{\partial{\bf{a}}} \\ \frac{\partial{v}}{\partial{\bf{a}}} \end{bmatrix} = \begin{bmatrix}\frac{\partial{u}}{\partial{w_x}} & \frac{\partial{u}}{\partial{w_y}} & \frac{\partial{u}}{\partial{w_z}} & \frac{\partial{u}}{\partial{f}} & \frac{\partial{u}}{\partial{u_0}} & \frac{\partial{u}}{\partial{v_0}}\\ \frac{\partial{v}}{\partial{w_x}} & \frac{\partial{v}}{\partial{w_y}} & \frac{\partial{v}}{\partial{w_z}} & \frac{\partial{v}}{\partial{f}} & \frac{\partial{v}}{\partial{u_0}} & \frac{\partial{v}}{\partial{v_0}}\end{bmatrix}$

In, general the expressions for these partial derivatives are quite complex. Most of the complexity arises from converting the rotation from axis-angle representation to rotation matrix representation. If we consider a 3D point that has already been transformed into the camera coordinate system, then the jacobian becomes significantly simpler and easier to analyze. Let’s denote the 3D point in the camera coordinate system by $P^\prime$ . $P^\prime$ is related to $P$ as:

$P^\prime = \begin{bmatrix}R & t\end{bmatrix}P$

The image coordinates $u$ and $v$ can be obtained by applying the perspective projection equation:

$u = \frac{fX^\prime}{Z^\prime}+u_0, v = \frac{fY^\prime}{Z^\prime}+v_0$

The jacobians are now easily obtained as:

$A = \begin{bmatrix}\frac{\partial{u}}{\partial{\bf{a}}} \\ \frac{\partial{v}}{\partial{\bf{a}}} \end{bmatrix} = \begin{bmatrix} \frac{\partial{u}}{\partial{f}} & \frac{\partial{u}}{\partial{u_0}} & \frac{\partial{u}}{\partial{v_0}}\\ \frac{\partial{v}}{\partial{f}} & \frac{\partial{v}}{\partial{u_0}} & \frac{\partial{v}}{\partial{v_0}}\end{bmatrix} = \begin{bmatrix} \frac{X^\prime}{Z^\prime} & 1 & 0 \\ \frac{Y^\prime}{Z^\prime} & 0 & 1\end{bmatrix}$

$B = \begin{bmatrix}\frac{\partial{u}}{\partial{\bf{b}}} \\ \frac{\partial{v}}{\partial{\bf{b}}} \end{bmatrix} = \begin{bmatrix} \frac{\partial{u}}{\partial{X^\prime}} & \frac{\partial{u}}{\partial{Y^\prime}} & \frac{\partial{u}}{\partial{Z^\prime}}\\ \frac{\partial{v}}{\partial{X^\prime}} & \frac{\partial{v}}{\partial{Y^\prime}} & \frac{\partial{v}}{\partial{Z^\prime}}\end{bmatrix} = \begin{bmatrix} \frac{f}{Z^\prime} & 0 & -\frac{fX}{Z^{\prime2}} \\ 0 & \frac{f}{Z^\prime} & -\frac{fY}{Z^{\prime2}}\end{bmatrix}$

Let’s look at the value of these jacobians at some reasonable values of $f, X^\prime, Y\prime, Z^\prime$ . Let’s consider a camera with a focal length of 5mm looking at a point located at $X^\prime=20, Y^\prime=20, Z^\prime=80$ . Using these values, the jacobian matrix is:

$B = \begin{bmatrix}\frac{\partial{u}}{\partial{\bf{b}}} \\ \frac{\partial{v}}{\partial{\bf{b}}} \end{bmatrix} = \begin{bmatrix} \frac{\partial{u}}{\partial{X^\prime}} & \frac{\partial{u}}{\partial{Y^\prime}} & \frac{\partial{u}}{\partial{Z^\prime}}\\ \frac{\partial{v}}{\partial{X^\prime}} & \frac{\partial{v}}{\partial{Y^\prime}} & \frac{\partial{v}}{\partial{Z^\prime}}\end{bmatrix} = \begin{bmatrix} 0.00625 & 0 & -0.0015625 \\ 0 & 0.00625 & -0.0015625\end{bmatrix}$

From this jacobian matrix, we can make the following observations.

The image coordinates are much more sensitive to the intrinsic parameters than to the coordinates of the object point in the camera coordinate system. Among the intrinsic parameters, the coordinates of the center of projection is the most important (for $\frac{X^\prime,Y^\prime}{Z^\prime}<1$ ). Thus, it is very important to calibrate the camera properly in practical applications of computer vision.
For $\frac{X^\prime,Y^\prime}{Z^\prime}<1$ , changes in $X^\prime$ have a larger effect than changes in $Z^\prime$ on the image coordinates. This agrees with our intuition that for small objects, movement along the viewing direction effects the image less than movement perpendicular to the viewing direction.

Recall that $\bf{P^\prime} = [X^\prime,Y^\prime, Z^\prime]$ are the coordinates of the object point in the camera coordinate system. Thus, $\bf{P^\prime} = [R | t]P$ . If the object and camera coordinates systems are aligned, $R = I$ and hence $[X^\prime,Y^\prime, Z^\prime] = [X, Y, Z] + [t_x, t_y, t_z]$ . For small objects, most of the contribution to $\bf{P^\prime}$ is from the camera-object translation vector. Thus, the effect of changes in the object point coordinates in the object coordinate system on the image coordinates is even smaller.

Let’s now consider the complete expressions for the Jacobians. Due to the complexity introduced by converting the axis-angle representation into a rotation matrix used in the perspective projection equation, calculating the jacobians by hand is very tedious. I use the symbolic toolbox in Matlab to compute the jacobians. Here’s the code:

syms u0 v0 f real
syms tx ty tz wx wy wz real
syms X Y Z real

% Ihe intrinsic parameter matrix
K=[f 0 u0;
    0 f v0;
    0 0 1];

% Expression for the rotation matrix based on the Rodrigues formula
theta=sqrt(wx^2+wy^2+wz^2);
omega=  [0 -wz wy;
         wz 0 -wx;
        -wy wx 0;];
R = eye(3) + (sin(theta)/theta)*omega + ((1-cos(theta))/theta^2)*(omega*omega);
% Expression for the translation vector
t=[tx;ty;tz];
% perspective projection of the model point (X,Y,Z)
uvs=K*[R t]*[X; Y; Z; 1];
u=uvs(1)/uvs(3);
v=uvs(2)/uvs(3);
% calculate the geometric distance in x and y direction
% u,v = the x and y positions of the projection of the corresponding model point
dx=u;
dy=v;
% Evaluate the symbolic expression of the Jacobian w.r.t. the estimated parameters
% Jx=jacobian(dx,[au,av,u0,v0,wx wy wz tx ty tz]);
% Jy=jacobian(dy,[au,av,u0,v0,wx wy wz tx ty tz]);
Jx_f = jacobian(dx,f);             
Jx_u0 = jacobian(dx,[u0]);        
Jx_v0 = jacobian(dx,[v0]);    
Jx_wx = jacobian(dx,[wx]);
Jx_wy = jacobian(dx,[wy]);
Jx_wz = jacobian(dx,[wz]);
Jx_tx = jacobian(dx,[tx]);
Jx_ty = jacobian(dx,[ty]);
Jx_tz = jacobian(dx,[tz]);
Jx_X = jacobian(dx,[X]);
Jx_Y = jacobian(dx,[Y]);
Jx_Z = jacobian(dx,[Z]);

Jy_f = jacobian(dy,f);
Jy_u0 = jacobian(dy,[u0]);
Jy_v0 = jacobian(dy,[v0]);
Jy_wx = jacobian(dy,[wx]);
Jy_wy = jacobian(dy,[wy]);
Jy_wz = jacobian(dy,[wz]);
Jy_tx = jacobian(dy,[tx]);
Jy_ty = jacobian(dy,[ty]);
Jy_tz = jacobian(dy,[tz]);
Jy_X = jacobian(dy,[X]);
Jy_Y = jacobian(dy,[Y]);
Jy_Z = jacobian(dy,[Z]);

% Convert to matlab functions that are stored as m files in the Jacobians directory. These functions
% can be called directly from other scripts.  
Jx_u0_func = matlabFunction(Jx_u0, 'File', 'Jacobians\Jx_u0');
Jx_v0_func = matlabFunction(Jx_v0, 'File', 'Jacobians\Jx_v0');
Jx_wx_func = matlabFunction(Jx_wx, 'File', 'Jacobians\Jx_wx');
Jx_wy_func = matlabFunction(Jx_wy, 'File', 'Jacobians\Jx_wy')
Jx_wz_func = matlabFunction(Jx_wz, 'File', 'Jacobians\Jx_wz')
Jx_tx_func = matlabFunction(Jx_tx, 'File', 'Jacobians\Jx_tx')
Jx_ty_func = matlabFunction(Jx_ty, 'File', 'Jacobians\Jx_ty')
Jx_tz_func = matlabFunction(Jx_tz, 'File', 'Jacobians\Jx_tz')
Jx_f_func = matlabFunction(Jx_f, 'File', 'Jacobians\Jx_f')
Jx_X_func = matlabFunction(Jx_X, 'File', 'Jacobians\Jx_X')
Jx_Y_func = matlabFunction(Jx_Y, 'File', 'Jacobians\Jx_Y')
Jx_Z_func = matlabFunction(Jx_Z, 'File', 'Jacobians\Jx_Z')

Jy_u0_func = matlabFunction(Jy_u0, 'File', 'Jacobians\Jy_u0');
Jy_v0_func = matlabFunction(Jy_v0, 'File', 'Jacobians\Jy_v0');
Jy_wx_func = matlabFunction(Jy_wx, 'File', 'Jacobians\Jy_wx');
Jy_wy_func = matlabFunction(Jy_wy, 'File', 'Jacobians\Jy_wy')
Jy_wz_func = matlabFunction(Jy_wz, 'File', 'Jacobians\Jy_wz')
Jy_tx_func = matlabFunction(Jy_tx, 'File', 'Jacobians\Jy_tx')
Jy_ty_func = matlabFunction(Jy_ty, 'File', 'Jacobians\Jy_ty')
Jy_tz_func = matlabFunction(Jy_tz, 'File', 'Jacobians\Jy_tz')
Jy_f_func = matlabFunction(Jy_f, 'File', 'Jacobians\Jy_f')
Jy_X_func = matlabFunction(Jy_X, 'File', 'Jacobians\Jy_X')
Jy_Y_func = matlabFunction(Jy_Y, 'File', 'Jacobians\Jy_Y')
Jy_Z_func = matlabFunction(Jy_Z, 'File', 'Jacobians\Jy_Z')

syms u0 v0 f real

syms tx ty tz wx wy wz real

syms X Y Z real

% Ihe intrinsic parameter matrix

K=[f 0 u0;

0 f v0;

0 0 1];

% Expression for the rotation matrix based on the Rodrigues formula

theta=sqrt(wx^2+wy^2+wz^2);

omega= [0 -wz wy;

wz 0 -wx;

-wy wx 0;];

R = eye(3) + (sin(theta)/theta)*omega + ((1-cos(theta))/theta^2)*(omega*omega);

% Expression for the translation vector

t=[tx;ty;tz];

% perspective projection of the model point (X,Y,Z)

uvs=K*[R t]*[X; Y; Z; 1];

u=uvs(1)/uvs(3);

v=uvs(2)/uvs(3);

% calculate the geometric distance in x and y direction

% u,v = the x and y positions of the projection of the corresponding model point

dx=u;

dy=v;

% Evaluate the symbolic expression of the Jacobian w.r.t. the estimated parameters

% Jx=jacobian(dx,[au,av,u0,v0,wx wy wz tx ty tz]);

% Jy=jacobian(dy,[au,av,u0,v0,wx wy wz tx ty tz]);

Jx_f = jacobian(dx,f);

Jx_u0 = jacobian(dx,[u0]);

Jx_v0 = jacobian(dx,[v0]);

Jx_wx = jacobian(dx,[wx]);

Jx_wy = jacobian(dx,[wy]);

Jx_wz = jacobian(dx,[wz]);

Jx_tx = jacobian(dx,[tx]);

Jx_ty = jacobian(dx,[ty]);

Jx_tz = jacobian(dx,[tz]);

Jx_X = jacobian(dx,[X]);

Jx_Y = jacobian(dx,[Y]);

Jx_Z = jacobian(dx,[Z]);

Jy_f = jacobian(dy,f);

Jy_u0 = jacobian(dy,[u0]);

Jy_v0 = jacobian(dy,[v0]);

Jy_wx = jacobian(dy,[wx]);

Jy_wy = jacobian(dy,[wy]);

Jy_wz = jacobian(dy,[wz]);

Jy_tx = jacobian(dy,[tx]);

Jy_ty = jacobian(dy,[ty]);

Jy_tz = jacobian(dy,[tz]);

Jy_X = jacobian(dy,[X]);

Jy_Y = jacobian(dy,[Y]);

Jy_Z = jacobian(dy,[Z]);

% Convert to matlab functions that are stored as m files in the Jacobians directory. These functions

% can be called directly from other scripts.

Jx_u0_func = matlabFunction(Jx_u0, 'File', 'Jacobians\Jx_u0');

Jx_v0_func = matlabFunction(Jx_v0, 'File', 'Jacobians\Jx_v0');

Jx_wx_func = matlabFunction(Jx_wx, 'File', 'Jacobians\Jx_wx');

Jx_wy_func = matlabFunction(Jx_wy, 'File', 'Jacobians\Jx_wy')

Jx_wz_func = matlabFunction(Jx_wz, 'File', 'Jacobians\Jx_wz')

Jx_tx_func = matlabFunction(Jx_tx, 'File', 'Jacobians\Jx_tx')

Jx_ty_func = matlabFunction(Jx_ty, 'File', 'Jacobians\Jx_ty')

Jx_tz_func = matlabFunction(Jx_tz, 'File', 'Jacobians\Jx_tz')

Jx_f_func = matlabFunction(Jx_f, 'File', 'Jacobians\Jx_f')

Jx_X_func = matlabFunction(Jx_X, 'File', 'Jacobians\Jx_X')

Jx_Y_func = matlabFunction(Jx_Y, 'File', 'Jacobians\Jx_Y')

Jx_Z_func = matlabFunction(Jx_Z, 'File', 'Jacobians\Jx_Z')

Jy_u0_func = matlabFunction(Jy_u0, 'File', 'Jacobians\Jy_u0');

Jy_v0_func = matlabFunction(Jy_v0, 'File', 'Jacobians\Jy_v0');

Jy_wx_func = matlabFunction(Jy_wx, 'File', 'Jacobians\Jy_wx');

Jy_wy_func = matlabFunction(Jy_wy, 'File', 'Jacobians\Jy_wy')

Jy_wz_func = matlabFunction(Jy_wz, 'File', 'Jacobians\Jy_wz')

Jy_tx_func = matlabFunction(Jy_tx, 'File', 'Jacobians\Jy_tx')

Jy_ty_func = matlabFunction(Jy_ty, 'File', 'Jacobians\Jy_ty')

Jy_tz_func = matlabFunction(Jy_tz, 'File', 'Jacobians\Jy_tz')

Jy_f_func = matlabFunction(Jy_f, 'File', 'Jacobians\Jy_f')

Jy_X_func = matlabFunction(Jy_X, 'File', 'Jacobians\Jy_X')

Jy_Y_func = matlabFunction(Jy_Y, 'File', 'Jacobians\Jy_Y')

Jy_Z_func = matlabFunction(Jy_Z, 'File', 'Jacobians\Jy_Z')

As mentioned earlier, the expressions for the jacobians are very complex, however most of the complexity arises from the axis-angle representation for rotation. For example, here’s the expression for Jy_z:

– (f*((wx*sin((wx^2 + wy^2 + wz^2)^(1/2)))/(wx^2 + wy^2 + wz^2)^(1/2) + (wy*wz*(cos((wx^2 + wy^2 + wz^2)^(1/2)) – 1))/(wx^2 + wy^2 + wz^2)) – v0*(((cos((wx^2 + wy^2 + wz^2)^(1/2)) – 1)*(wx^2 + wy^2))/(wx^2 + wy^2 + wz^2) + 1))/(tz – X*((wy*sin((wx^2 + wy^2 + wz^2)^(1/2)))/(wx^2 + wy^2 + wz^2)^(1/2) + (wx*wz*(cos((wx^2 + wy^2 + wz^2)^(1/2)) – 1))/(wx^2 + wy^2 + wz^2)) + Y*((wx*sin((wx^2 + wy^2 + wz^2)^(1/2)))/(wx^2 + wy^2 + wz^2)^(1/2) – (wy*wz*(cos((wx^2 + wy^2 + wz^2)^(1/2)) – 1))/(wx^2 + wy^2 + wz^2)) + Z*(((cos((wx^2 + wy^2 + wz^2)^(1/2)) – 1)*(wx^2 + wy^2))/(wx^2 + wy^2 + wz^2) + 1)) – ((((cos((wx^2 + wy^2 + wz^2)^(1/2)) – 1)*(wx^2 + wy^2))/(wx^2 + wy^2 + wz^2) + 1)*(f*ty + X*(f*((wz*sin((wx^2 + wy^2 + wz^2)^(1/2)))/(wx^2 + wy^2 + wz^2)^(1/2) – (wx*wy*(cos((wx^2 + wy^2 + wz^2)^(1/2)) – 1))/(wx^2 + wy^2 + wz^2)) – v0*((wy*sin((wx^2 + wy^2 + wz^2)^(1/2)))/(wx^2 + wy^2 + wz^2)^(1/2) + (wx*wz*(cos((wx^2 + wy^2 + wz^2)^(1/2)) – 1))/(wx^2 + wy^2 + wz^2))) + tz*v0 + Y*(v0*((wx*sin((wx^2 + wy^2 + wz^2)^(1/2)))/(wx^2 + wy^2 + wz^2)^(1/2) – (wy*wz*(cos((wx^2 + wy^2 + wz^2)^(1/2)) – 1))/(wx^2 + wy^2 + wz^2)) + f*(((cos((wx^2 + wy^2 + wz^2)^(1/2)) – 1)*(wx^2 + wz^2))/(wx^2 + wy^2 + wz^2) + 1)) – Z*(f*((wx*sin((wx^2 + wy^2 + wz^2)^(1/2)))/(wx^2 + wy^2 + wz^2)^(1/2) + (wy*wz*(cos((wx^2 + wy^2 + wz^2)^(1/2)) – 1))/(wx^2 + wy^2 + wz^2)) – v0*(((cos((wx^2 + wy^2 + wz^2)^(1/2)) – 1)*(wx^2 + wy^2))/(wx^2 + wy^2 + wz^2) + 1))))/(tz – X*((wy*sin((wx^2 + wy^2 + wz^2)^(1/2)))/(wx^2 + wy^2 + wz^2)^(1/2) + (wx*wz*(cos((wx^2 + wy^2 + wz^2)^(1/2)) – 1))/(wx^2 + wy^2 + wz^2)) + Y*((wx*sin((wx^2 + wy^2 + wz^2)^(1/2)))/(wx^2 + wy^2 + wz^2)^(1/2) – (wy*wz*(cos((wx^2 + wy^2 + wz^2)^(1/2)) – 1))/(wx^2 + wy^2 + wz^2)) + Z*(((cos((wx^2 + wy^2 + wz^2)^(1/2)) – 1)*(wx^2 + wy^2))/(wx^2 + wy^2 + wz^2) + 1))^2

🙂

Lourakis MIA, Argyros AA. SBA. ACM Transactions on Mathematical Software. 2009;36(1):1-30. doi: 10.1145/1486525.1486527 [Source]

Lourakis MIA, Argyros AA. SBA. ACM Transactions on Mathematical Software. 2009;36(1):1-30. doi:10.1145/1486525.1486527

Dao Khanh Hoai
September 18, 2017 at 4:20 am

Hello Ankur!
Your works has really impressed me. I’m a researcher from Ha Noi – Viet Nam. I have read and tryed to replicate your steps in Bundle Adjustment. I would like you to do a favor: Share your simulated data sets? Or You show how to create these data sets.
Thank you very much!

lyubomir dinchev
March 1, 2018 at 10:46 am

Grate article, you can use first order approximation for the exponential mapping of the Lie group/algebra, it simplifies quite a lot the Jacobian. Here is an example: https://vision.in.tum.de/_media/spezial/bib/kerl2012msc.pdf

- lyubomir dinchev
  March 1, 2018 at 10:51 am
  
  And I can learn to not write “Grate” when i want to write “Great”.

Telesens

Bundle Adjustment – Part 1: Jacobians

3 Comments

Leave a Reply Cancel reply