Need help for proper rectilinear projection of equirectangular panoramic image

Question

With the algorithm below, when the projection plane is tangent to the equator (the center line of the equirectangular image), projected image looks rectilinear.

But when the projection plane is tilted, (py0 != panorama.height/2), lines are warped.

The two last "lines" in the algorithm below needs to be "rectified", in order to adjust px and/or py when the center line of the destination plane is not at the same level than the center line of the equirectangular image.

  // u,v,w :
  //     Normalized 3D coordinates of the destination pixel

  // elevation, azimuth:
  //     Angles between the origin (sphere center) and the destination pixel

  // px0, py0 :
  //     2D coordinates in the equirectangular image for the
  //     the destination plane center (long*scale,lat*scale)

  // px, py:
  //     2D coordinates of the source pixel in the equirectangular image
  //     (long*scale,lat*scale)

  angularStep=2*PI/panorama.width;
  elevation=asin(v/sqrt(u*u+v*v+w*w));
  azimuth=-PI/2+atan2(w,u);
  px=px0+azimuth/angularStep;
  py=py0+elevation/angularStep;

I can compute the intersection p between the normal of each destination pixel and the sphere, then convert cartesian coordinates to long/lat using available C code:

But I know there's a simpler and much less time consuming method, involving adjusting source pixel coordinates in equirectangular image (px,py) knowing the longitude/latitude (px0,py0) at which the center of the projection plane intersect the "sphere".

Could you help please ?

Answer 1

I managed to get this to work using the formula for gnomonic projection in a webgl shader http://mathworld.wolfram.com/GnomonicProjection.html

            float angleOfView   

            float phi1          
            float lambda0     //centre of output projection

            float x = PI2*(vTextureCoord.s - 0.5) ;  //input texture coordinates, 
            float y = PI2*(vTextureCoord.t - 0.5 ); 

            float p = sqrt(x*x + y*y);

            float c = atan(p, angleOfView); 

            float phi = asin( cos(c)*sin(phi1) + y*sin(c)*cos(phi1)/p );

            float lambda = lambda0 + atan( x*sin(c), (p*cos(phi1)*cos(c) - y*sin(phi1)*sin(c)));

            vec2 tc = vec2((lambda /(PI*2.0) + 0.5, (phi/PI) + 0.5); //reprojected texture coordinates 

            vec4 texSample  =  texture2D(tEqui, tc); //sample using new coordinates