How to find locations of 3D points projected to panorama? (rotation order problem)

Question

I'm trying to find the 2D position of where a 3D point will be when projected to an equirectangular panoramic render.

I have a function for doing this with perspective cameras but this doesn't work when the camera type is panorama (even in cycles), probably because panoramas don't have a linear projection matrix:

def project(cam, point):
    co_2d = bpy_extras.object_utils.world_to_camera_view(C.scene, cam, point)
    render_scale = C.scene.render.resolution_percentage / 100
    render_size = (int(C.scene.render.resolution_x * render_scale), int(C.scene.render.resolution_y * render_scale))
    return (int(co_2d.x * render_size[0]) % C.scene.render.resolution_x, int(co_2d.y * render_size[1]) % C.scene.render.resolution_y)

So I went about trying to write the function from scratch and got most of the way.. In an equirectangular image the x and y axes correspond to longitude and latitude respectively. Thus to find the screen position we just need to convert the vector between the camera and target point into spherical polars:

def polar_global(camera, point):
    v_p = Vector((camera.location - point)).normalized()
    x, y, z = v_p[0:3]
    longitude = (pi - atan2(y, x)) * renX/(2*pi)
    latitude = ((pi/2) - atan2(z, sqrt(x*x + y*y))) * renY/pi
    return (int(longitude), int(latitude))

This works correctly for any camera location, i.e. the calculated positions line up with a target panoramic render where the camera is pointing down the x axis, i.e. camera.rotation_euler = (90,0,90), for any camera location:

Now the problem is when I try to incorporate the cameras rotation. This should be as simple as rotating each vector by the inverse of the cameras rotation before converting to polar coordinates, however because the camera is offset by (90,0,90) to be pointing along the x axis, this must be corrected before multiplying to get the relative vector.

I've tried a few combinations such as negating the components of the rotation, etc. but I cannot get it to line up. I've been subtracting radians(90) from the x and z components before inverting the rotation since this seemed to make more sense but it might need to be done after.

Also doing say, z - 90 instead of 90 - z also gives a rotation of 0 when z is 90 (pointing down the axis) so its confusing which to use.

def polar(camera, point):
v_p = Vector((camera.location - point)).normalized()
cam_E = Euler((-camera.rotation_euler[0] - radians(90), -camera.rotation_euler[1], -camera.rotation_euler[2] - radians(90)), 'XYZ') 
v_p_r = cam_E.to_quaternion().inverted() @ v_p
x, y, z = v_p_r[0:3]
longitude = (pi - atan2(y, x)) * renX/(2*pi)
latitude = ((pi/2) - atan2(z, sqrt(x*x + y*y))) * renY/pi
return (int(longitude), int(latitude))

Blend file: blend

Answer 1

After a lot of thinking (and a hot bath), I figured I couldn't just subtract the rotation offset in the Euler rotation so I separated out this out and treated them as separate rotations. I tried a few different orders and the one that worked agreed with what you'd expect using matrix math..

The camera rotation in the scene is made up of the offset rotation first, followed by the rotation we desire:

camera = target @ offset

Hence to get the desired camera rotation the offset must be undone before the camera rotation: target = camera @ offset.inverted()

Finally, this target is in fact how much the virtual camera has rotated, and so we must invert the whole thing so that the vectors to our points in the scene rotate the opposite way to this:

rotation = (camera @ offset.inverted()).inverted() = offset @ camera.inverted()

Final code:

def polar(camera, point):
    v_p = Vector((camera.location - point)).normalized()
    d_cam = camera.rotation_euler.to_quaternion()
    d_offset = Euler((radians(90), 0, radians(90)), 'XYZ').to_quaternion()
    v_p_r = d_offset @ d_cam.inverted() @ v_p
    x, y, z = v_p_r[0:3]
    longitude = (pi - atan2(y, x)) * renX/(2*pi)
    latitude = ((pi/2) - atan2(z, sqrt(x*x + y*y))) * renY/pi
    return (int(longitude), int(latitude))