Using Image Transformation to create image sky view of what its like to be near black hole (Not a wormhole)

Question

Question : I am trying to carry out point 3. "Implementing map" on page 13 in the paper Interstellar wormholes using my own equations not those in the paper. however am not 100% sure on how to go about generating the numerical table of values and the image processing parts.

Note: I am using different equations to those given in the paper and thus the result will be different however the general method should be similiar therefore this is NOT the same answer as

I have generated the numerical map all I am trying to do is put this into the ImageTransformation or ImageFowardTransformation function, not entirely sure how to relate my interpolation function to the angles of a 2d image

numericalmap = {};
n = 100;
For[i = 0, i < n + 1, i++, \[Phi]csgen = (1.009 + (1/2) i/n)*Pi;
 M = 1; E0 = 1; \[Theta]cs = Pi/2;
 Nx = Sin[\[Theta]cs]*Cos[\[Phi]csgen];
 Ny = Sin[\[Theta]cs]*Sin[\[Phi]csgen];
 Nz = Cos[\[Theta]cs];
 nr = -Nx;
 n\[Phi] = -Ny;
 n\[Theta] = Nz;
 b = rc*Sin[\[Theta]c]*n\[Phi]/(1 - 2 M/rc)^(1/2);
 B2 = rc^2/(1 - 2 M/rc)*(n\[Phi]^2 + n\[Theta]^2);
 prinitial = ((1 - 2 M/rc)^(-1))*nr;
 p\[Theta]initial = ((1 - 2 M/rc)^(-1/2))*rc*n\[Theta];
 {rc, \[Theta]c, \[Phi]c} = {200, Pi/2, 0};
 lambdaend = -100000;
 ham = {
   t'[\[Lambda]] + E0/(1 - (2 M)/r[\[Lambda]]) == 0,
   r'[\[Lambda]] - (1 - (2 M)/r[\[Lambda]]) pr[\[Lambda]] == 0,
   \[Theta]'[\[Lambda]] - P\[Theta][\[Lambda]]/r[\[Lambda]]^2 == 0,
   \[Phi]'[\[Lambda]] - b/(r[\[Lambda]]*Sin[\[Theta][\[Lambda]]])^2 ==
     0, pr'[\[Lambda]] + 
     M/r[\[Lambda]]^2 (E0^2/(1 - (2 M)/r[\[Lambda]])^2 + 
        pr[\[Lambda]]^2) - B2/r[\[Lambda]]^3 == 0, 
   P\[Theta]'[\[Lambda]] - (b^2*
        Cos[\[Theta][\[Lambda]]])/(r[\[Lambda]]^2*
        Sin[\[Theta][\[Lambda]]]^3) == 0
   };
 haminital = {t[0] == 0, 
   r[0] == rc, \[Theta][0] == \[Theta]c, \[Phi][0] == \[Phi]c, 
   pr[0] == prinitial, P\[Theta][0] == p\[Theta]initial};
 \[Phi]2 = 
  NDSolveValue[{ham, haminital}, {t, r, \[Theta], \[Phi], pr, 
     P\[Theta]}, {\[Lambda], lambdaend, 0}][[4]];
 numericalmap = 
  Append[numericalmap, {\[Phi]csgen - Pi, Pi + \[Phi]2[lambdaend]}]]

Would appreciate any insight as how I could input this interpolation into ImageTransformation to get the desired effect.

I know Jason B attempted a similar problem but didn't use ImageTransformation for it.

When I put in some image of space the resulting images should look something like this: