Here is a proof for the plane. I think the same method of proof directly generalizes to higher dimensions, but it might get annoying to write down.
DEFINITION: A projectivity is a function \varphi on \mathbb{R}^2 such that
\varphi(x,y)=\left(\frac{Ax+By+C}{Gx+Hy+I},\frac{Dx+Ey+F}{Gx+Hy+I}\right)
where A,B,C,D,E,F,G,H,I are real numbers such that the matrix
\left(\begin{array}{ccc} A & B & C\\ D & E & F\\ G & H & I\end{array}\right)
is invertible. This invertible-condition tells us exactly that \varphi is invertible. The inverse is again a perspectivity and its matrix is given by the inverse of the above matrix.
We can see this easily as follows:
Recall that a homogeneous coordinate is defined as a triple [x:y:z] with not all x, y and z zero. Furthermore, if \alpha\neq 0, then we define [\alpha x: \alpha y : \alpha z]=[x:y:z].
There exists a bijection between \mathbb{R}^2 and the homogeneous coordinates [x:y:z] with nonzero z. Indeed, with (x,y) in \mathbb{R}^2, we can associate [x:y:1]. And with [x:y:z] with nonzero z, we can associate (x/z,y/z).
We can now look at \varphi on homogeneous coordinates. We define \varphi [x:y:z] = \varphi(x/z,y/z). Clearly, if \alpha\neq 0, then \varphi [\alpha x:\alpha y:\alpha z]=\varphi [x:y:z]. So the map is well defined.
Actually, our \varphi is actually just matrix multiplication:
\varphi[x:y:z] = \left(\begin{array}{ccc} A & B & C\\ D & E & F\\ G & H & I\end{array}\right)\left(\begin{array}{c} x\\ y \\ z\end{array}\right)
Now we see clearly that \varphi has an inverse given by
\varphi^{-1} [x:y:z] = \left(\begin{array}{ccc} A & B & C\\ D & E & F\\ G & H & I\end{array}\right)^{-1}\left(\begin{array}{c} x\\ y \\ z\end{array}\right)
LEMMA: Let x,y,z and t in \mathbb{R}^2 be four distinct points such that no three of them lie on the same line. Let x',y',z',t' in \mathbb{R}^2 also be four points such that no three of them lie on the same line. There exists a projectivity \varphi such that \varphi(x)=x^\prime, \varphi(y)=y^\prime, \varphi(z)=z^\prime, \varphi(t)=t^\prime.
We write in homogeneous coordinates:
x=[x_1:x_2:x_3],~y=[y_1:y_2:y_3],~z=[z_1:z_2:z_3],~t=[t_1:t_2:t_3]
Since \mathbb{R}^3 has dimension 3, we can find \alpha,\beta,\gamma in \mathbb{R} such that
(t_1,t_2,t_3)=(\alpha x_1,\alpha x_2,\alpha x_3)+(\beta y_1,\beta y_2,\beta y_3)+ (\gamma z_1, \gamma z_2,\gamma z_3).
The vectors (\alpha x_1,\alpha x_2,\alpha x_3), (\beta y_1,\beta y_2,\beta y_3), (\gamma z_1, \gamma z_2,\gamma z_3) form a basis for \mathbb{R}^3 (because of the condition that not three of x,y,z or t is on one line).
We can do the same for the x',y',z',t' and we again obtain a basis (\alpha^\prime x_1^\prime,\alpha^\prime x_2^\prime,\alpha^\prime x_3^\prime), (\beta^\prime y_1^\prime,\beta^\prime y_2^\prime,\beta^\prime y_3^\prime), (\gamma^\prime z_1^\prime, \gamma^\prime z_2^\prime,\gamma^\prime z_3^\prime) such that
(t_1^\prime, t_2^\prime,t_3^\prime)=(\alpha^\prime x_1^\prime,\alpha^\prime x_2^\prime,\alpha^\prime x_3^\prime)+(\beta^\prime y_1^\prime,\beta^\prime y_2^\prime,\beta^\prime y_3^\prime)+(\gamma^\prime z_1^\prime, \gamma^\prime z_2^\prime,\gamma^\prime z_3^\prime)
By linear algebra, we know that there exists an invertible matrix T that sends the bases on each other. This implies directly that the associated projectivity sends x to x', y to y' and z to z'.
Since
(t_1,t_2,t_3)=(\alpha x_1,\alpha x_2,\alpha x_3)+(\beta y_1,\beta y_2,\beta y_3)+ (\gamma z_1, \gamma z_2,\gamma z_3)
we get after applying T that
T(t_1,t_2,t_3)=(\alpha^\prime x_1^\prime,\alpha^\prime x_2^\prime,\alpha^\prime x_3^\prime)+(\beta^\prime y_1^\prime,\beta^\prime y_2^\prime,\beta^\prime y_3^\prime)+(\gamma^\prime z_1^\prime, \gamma^\prime z_2^\prime,\gamma^\prime z_3^\prime)
and thus T(t_1,t_2,t_3)=(t_1^\prime,t_2^\prime, t_3^\prime). Thus the projectivity also sends t to t'.
THEOREM Let U\subseteq \mathbb{R}^2 be open and let \varphi:U\rightarrow \mathbb{R}^2 be injective. Assume that \varphi sends lines to lines, then it is a projectivity.
We can of course assume that U contains an equilateral triangle ABC. Let P be the centroid of ABC.
By the previous lemma, there exists a projectivity \psi such that \psi(\varphi(A))=A, ~\psi(\varphi(B))=B, ~\psi(\varphi(C))=C, ~\psi(\varphi(P))=P. So we see that \sigma:=\psi\circ\varphi sends lines to lines and that \sigma(A)=A,~\sigma(B)=B,~\sigma(C)=C,~\sigma(P)=P. We will prove that \sigma is the identity.
HINT: look at Figure 2.1, p.19 of the Mccallum paper.
Define E the midpoint of AC. Then E is the intersection of AC and PB. But these lines are fixed by \sigma. Thus \sigma(E)=E. Let D be the midpoint of BC and F the midpoint of AB. Likewise follows that \sigma(D)=D and \sigma(F)=F.
Thus \sigma preserves the verticles of the equilateral traingles AFE, FBD, DEF and EDC. Since \sigma preserves parallelism, we see easily that \sigma preserves the midpoints and centroids of the smaller triangles. So we can subdivide the triangles in even smaller triangles whose vertices are preserved. We keep doing this process and eventually we find a set S dense in the triangle such that \sigma is fixed on that dense set. If \sigma were continuous, then \sigma is the identity on the triangle.
To prove continuity, we show that certain rhombuses are preserved. Look at Figure 2.3 on page 20 of McCallum. We have shown that the vertices of arbitrary triangles are preserved. Putting those two triangles together gives a rhombus. We will show that \sigma sends the interior of any rhombus ABCD into the rhombus ABCD. Since the rhombus can be made arbitrarily small around an arbitrary point, it would follow that \sigma were continuous.
By composing with a suitable linear map, we restrict to the following situation:
LEMMA: Let A=(0,0), B=(1,0), C=(1,1) and D=(0,1) and let \Sigma be the square ABCD. Suppose that \sigma:\Sigma\rightarrow \mathbb{R}^2 sends lines to lines and suppose that \sigma is fixed on A,B,C and D. Then \sigma(\Sigma)\subseteq \Sigma.
Take S on CB. We can make a construction analogous to 2.4 p.22 in MCCallen. So we let TS be horizontal, TU have slope -1 and VU be vertical. We define Q as the intersection of AS and VU. If S has coordinates (1,s) for some s. Then we can easily check that Q has coordinates (s,s^2). In particular, Q lies in the upper half-plane (= everything about AB).
Since S in CB and since C and B are fixed. We see that \sigma(S)\in CB. Let's say that \sigma(S)=(1,t) for some t. The line TS is a horizontal and \sigma maps this to a horizontal. So \sigma(T) has the form (0,t). The line TU has slope -1. So \sigma(U) has the form (t,0). Finally, it follows that \sigma(Q) has the form (t,t^2). In particular, \sigma(Q) is in the upper half plane.
So we have proven that if S is on CB, then they ray AS emanating from A is sent into the upper half plane. Let P be an arbitrary point in the square, then it is an element of a ray AS for some S. This ray is taken to the upper half plane. So \sigma(P) is in the upper half plane.
So this proves that the square ABCD is sent by \sigma into the upper half plane. Similar constructions show that the square is also sent to the lower half plane, the left and right half planes. So taking all of these things together: ABCD is sent into ABCD. This proves the lemma.
So, right now we have shown that \sigma is the identity on some small equilateral triangle in U. So \varphi is a projectivity on some small open set U^\prime of U (namely on the interior of the triangle). We prove now that \varphi will be a projectivity on entire U.
Around any point P in U, we can find some equilateral triangle. And we proved for such triangles that \varphi is a projectivity and thus analytic. The uniqueness of analytic continuation now proves that \varphi is a projectivity on entire U.