I'm digitising some hand-drawn maps. My method is to scan the map, import it into an SVG editor and trace over the features. I then mark some datum points on the SVG, with metadata to indicate that they are datum points. Then I find those same datum points on Google Maps to get their location in decimal lat/long and add those values to the metadata. Note that the hand-drawn map is not necessarily oriented so that the Y axis points North.
So now I have two lists of vectors: $A$ is an 2xN matrix with the datum coordinates in SVG coordinates as $(x, y)$ and $B$ is a 2xN matrix with the datum coordinates as $(lng, lat)$.
I now want to find a rotation, scale and translation that will map coordinates from one space to the other, so that I can then transform the actual features on the map to lat/long.
I found this method which uses SVD to obtain a rotation and translation matrix, but I'm not sure how well the SVD step will work when there is a scale involved as well. To summarise the method:
- Let $A_c$ and $B_c$ be the centroids of $A$ and $B$ (ie the arithmetic mean along each column).
- Calculate the covariance matrix $H=(A-A_c)(B-B_c)^T$.
- Do the SVD -- $U, S, V = SVD(H)$.
- The rotation matrix is then $R = VU^T$.
- The translation is $t=B-RA$
Can I then just calculate the scale as the mean along each column $S = \overline{B_c / (RA_c)}$ (where / is element-wise division)? And then re-evaluate the translation as $t=B-diag(S)(RA)$? Or does the scale need to be taken into account in the SVD somehow?
The same website that you linked has an article that explains the inner workings of OpenCV's estimateRigidTransform function. I believe that's what you are looking for to compute the scale.
http://nghiaho.com/?p=2208
Assuming uniform scale (no shear), solve for [a, b, c, d] using 2 points: $$\begin{bmatrix} x_1 & -y_1 & 1 & 0 \\ y_1 & x_1 & 0 & 1 \\ x_2 & -y_2 & 1 & 0 \\ y_2 & x_2 & 0 & 1 \end{bmatrix} \begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix} = \begin{bmatrix} x’_1 \\ y’_1 \\ x’_2 \\ y’_2 \end{bmatrix}$$
Where: $$a = \cos(\theta)s$$ $$b = \sin(\theta)s$$ $$c = tx$$ $$d = ty$$
If you have 3 or more points, a simple least square solution can be obtained by doing a pseudo-inverse: $$\begin{matrix} Ax & = & b \\ A{^T}Ax & = & A^{T}b \\x & =& (A{^T}A)^{-1}A{^T}b \end{matrix}$$