Plane through centre of unit sphere appears offset in OpenGL ES 2.0

Question

I'm plotting some coastline data on a sphere using the vispy interface for OpenGL ES 2.0. I'm using latitude and longitude values to work out 3d coordinates for the data on a sphere and just plotting these. I'm able to successfully draw the data but I wanted to see only those data points on the side of the sphere that would be visible from the viewport.

I've tried two quite different approaches to create this effect but both have led to the same problem. First, I calculated the dot product of the view direction and data position and drew only those with a negative result (i.e. only those points facing the viewport) and, secondly, I simply drew a plane through the centre of the sphere, perpendicular to the view direction.

In both cases I observed the same - the plane appeared to be slightly offset away from the viewport, behind the centre of the sphere. In other words, you can see the data wrap around the back of the sphere slightly before it's masked by the plane.

I've checked that the points I'm drawing are, in fact, on the unit sphere and I feel confident that, from the 3d world point of view, everything is sound. What I am less confident with, as a relative beginner to 3d graphics, is whether I'm misunderstanding something with the projection matrix. I've done some reading - but my understanding leads me to think that the projection shouldn't change the order of points in the "Z direction" (the direction the viewport is facing).

I'm confident this isn't a depth test issue as my first approach didn't have depth test enabled and masking was done in the vertex shader (by setting the fragment colour alpha to 0.0). Aside from this, I've not been able to find any other explanation for the issue.

Here's the code for the plane approach:

import numpy as np
import cartopy
from vispy import app
from vispy import gloo
import time
from vispy.util.transforms import perspective, translate, rotate

xpts = []
ypts = []

#getting coastlines data

for string in cartopy.feature.NaturalEarthFeature('physical', 'coastline', '10m').geometries():
    for line in string:
        points = list(line.coords)
        for point in points:
            xpts.append(point[0])
            ypts.append(point[1])

coasts = np.array(zip(xpts,ypts), dtype=np.float32)

theta = (np.pi/180)*np.array(xpts, dtype=np.float32)
phi = (np.pi/180)*np.array(ypts, dtype=np.float32)

x3d = np.cos(phi)*np.cos(theta)
y3d = np.sin(theta)*np.cos(phi)
z3d = np.sin(phi)



vertex = """
// Uniforms
uniform mat4 u_model;
uniform mat4 u_view;
uniform mat4 u_projection;
uniform vec3 u_color;

attribute vec3 a_position;
void main (void)
{
    gl_Position = u_projection*u_view*u_model*vec4(a_position, 1.0);
}
"""

fragment = """
// Uniforms
uniform vec3 u_color;

void main()
{
    gl_FragColor = vec4(u_color, 1.0);
}
"""

class Canvas(app.Canvas):
    def __init__(self):
        app.Canvas.__init__(self, keys='interactive')

        gloo.set_state(clear_color = 'red', depth_test=True, blend=True, blend_func=('src_alpha', 'one_minus_src_alpha'))

        self.x = 0

        self.plane = 5*np.array([(0.,-1., -1.,1), (0, -1., +1.,1), (0, +1., -1.,1), (0, +1., +1.,1)], dtype=np.float32)

        self._timer = app.Timer(connect=self.on_timer, start=True)
        self.program = gloo.Program(vertex, fragment)

        self.view = np.dot(rotate(-90, (1, 0, 0)), np.dot(translate((-3, 0, 0)), rotate(-90.0, (0.0,1.0,0.0))))
        self.model = np.eye(4, dtype=np.float32)
        self.projection = perspective(45.0, self.size[0]/float(self.size[1]), 2.0, 10.0)

        self.program['u_projection'] = self.projection
        self.program['u_view'] = self.view
        self.program['u_model'] = self.model
        self.program['u_color'] = np.array([0.0, 0.0, 0.0], dtype=np.float32)

        self.program2 = gloo.Program(vertex, fragment)

        self.program2['u_projection'] = self.projection
        self.program2['u_view'] = self.view
        self.program2['u_model'] = self.model
        self.program2['u_color'] = np.array([1.0, 1.0, 1.0], dtype=np.float32)

        self.program2['a_position'] =  self.plane[:,:3].astype(np.float32)

    def on_timer(self, event):
        self.x += 0.05
        self.model = rotate(self.x, (0.0,0.0,1.0))
        pointys = np.concatenate((x3d,y3d,z3d)).reshape((3, -1)).T
        self.program['a_position'] = pointys
        self.program['u_model'] = self.model

        self.update()


    def on_resize(self, event):
        gloo.set_viewport(0, 0, *event.size)
        self.projection = perspective(45.0, event.size[0]/float(event.size[1]), 2.0, 10.0)
        self.program['u_projection'] = self.projection
        self.program2['u_projection'] = self.projection


    def on_draw(self, event):
        gloo.clear((1,1,1,1))
        self.program2.draw('triangle_strip')
        self.program.draw('points')

Canvas().show()
app.run()

Answer 1

The way I understand your description, what you're seeing is a result of the perspective projection. I used all of my MS Paint skills to create this very elaborate diagram of the situation viewed from the side:

The outline of the sphere is drawn in black. The red line indicates a plane through the center of the sphere.

The blue lines show two lines of sight from the viewpoint, which is at the bottom of the diagram. If you picture the result after applying the projection, what shows up as the front facing part of the sphere in the rendered image is everything below the green line. The parts of the sphere above the green line form the back facing part of the sphere in the resulting rendering.

Or in other words, the green line shows the plane that corresponds to the outline of the sphere in the resulting rendering.

As you can see from this, the plane through the center of the sphere is indeed some distance behind the section of the sphere that shows up as the front facing part of the sphere in the rendered image. This is just in the nature of a perspective projection. The distance between the red plane and the green plane will decrease with a smaller viewing angle (i.e. a weaker perspective), and the two are the same when using a parallel projection.