Let's say the coordinates of the pixels are $(\theta_i, \phi_i)$. $\theta$ is measured from the equator, so it varies between $\pm \pi/2$ and $\phi$ is measured in the horizontal plane, from $0$ to $2\pi$. This is like the standard spherical coordinate system shifted by $\pi/2$. Then the unit vectors in Cartesian coordinates for a point are described by:
$$\begin{align}x_i&=\cos\theta_i\cos\phi_i\\y_i&=\cos\theta_i\sin\phi_i\\z_i&=\sin\theta_i\end{align}$$
As mentioned in another answer, we need to find a rotation axis. For that, we use the cross product:
$$\vec n=\vec P_1\times\vec P_2$$
By components, this becomes: $$\begin{align}n_x&=y_1z_2-z_1y_2\\n_y&=z_1x_2-x_1z_2\\n_z&=x_1y_2-y_1x_2\end{align}$$
Note that $\vec n$ doess not necessarily have length $1$, so you will need to use $\hat n=\frac{\vec n}{|\vec n|}$. Also, $\vec n$ might be $0$. Then you have two points either identical or on diametrically opposite positions. In that case there are infinitely many rotation axes you can choose, in the perpendicular plane to the diameter. So for example, if the two points are the North and South poles, then any meridian would work. In that case, I would calculate the angle between the diameter and each of the $x, y, z$ axes, choose the largest one of them (or the first of the largest, if they are equal), and use the Gram-Schmidt orthogonalization process to find $\hat n$. The reason I would choose the largest is that it will introduce less numerical errors.
Note that $|\vec n|=|\vec P_1||\vec P_2|\sin\alpha=\sin\alpha$., where $\alpha$ is the angle between the vectors. If $|\vec n|=0$, you need to check if the points are identical $\alpha =0$ or opposite $\alpha =\pi$.
With the given information, we can now proceed to find the intermediate points. for that, we create an array of $\alpha_i$ from $0$ to $\alpha$, and we rotate the vector $\vec P_1$ towards vector $\vec P_2$ in the plane perpendicular to $\hat n$ using all $\alpha_i$ rotation angles. You can use the formula on wikipedia for the rotation matrix:
$$R(\hat n, \alpha_i)=\begin{pmatrix}\cos\alpha_i+u_x^2(1-\cos\alpha_i)&u_xu_y(1-\cos\alpha_i)-u_z\sin\alpha_i&u_xu_z(1-\cos\alpha_i)+u_y\sin\alpha_i\\u_yu_x(1-\cos\alpha_i)+u_z\sin\alpha_i&\cos\alpha_i+u_y^2(1-\cos\alpha_i)&u_yu_z(1-\cos\alpha_i)-u_x\sin\alpha_i\\u_xu_z(1-\cos\alpha_i)-u_y\sin\alpha_i&u_zu_y(1-\cos\alpha_i)+u_x\sin\alpha_i&\cos\alpha_i+u_z^2(1-\cos\alpha_i)\end{pmatrix}$$
I've used the components of $\hat n$: $$u_j=\frac{n_j}{|\vec n|}$$Then the intermediate points are $$\vec P_i=R(\hat n,\alpha_i)\vec P_1$$This is a simple matrix multiplication for each point. You can now calculate back the $\theta_i$ and $\phi_i$ coordinates, if you want to use those.
