Order of vertex for drawing quad

Question

I need to draw n planes in 3D space. The quads are planes created from two points, and with an algorithm I get 4 vertex for drawing the quad. The problem I have is that the order of the vertices affect to the result, obviously. Here it is what I mean: But when the plane is horizontal instead of vertical:

I can imagine two possible solutions. Using triangles and combining them or ordering the vertex properly. I dont know how to do the second idea. I have tried using triangles but I get the same problem.

# self.planos = [('A', (500, 500, 10), (-500, 500, 10), (-500, -500, 10), (500, -500, 10))] for horizontal
# self.planos = [('A', (-500, 10, 500), (500, 10, 500), (-500, 10, -500), (500, 10, -500))] for vertical
glEnable(GL_BLEND)
glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)
glDepthMask(GL_FALSE)
glBegin(GL_QUADS)
glColor(0.5, 0.5, 0.1, 0.5)
for i in range(len(self.planos)):
    glVertex(self.planos[i][1][0], self.planos[i][1][2], self.planos[i][1][1])
    glVertex(self.planos[i][2][0], self.planos[i][2][2], self.planos[i][2][1])
    glVertex(self.planos[i][3][0], self.planos[i][3][2], self.planos[i][3][1])
    glVertex(self.planos[i][4][0], self.planos[i][4][2], self.planos[i][4][1])
glEnd()
glDepthMask(GL_TRUE)
glDisable(GL_BLEND)

Intersection code for getting four vertex to draw planes:

In init method:
#Vertices of the cube
self.v = (Point3D(500, 500, 500), Point3D(-500, 500, 500), Point3D(-500, -500, 500),
          Point3D(500, -500, 500), Point3D(500, 500, -500), Point3D(-500, 500, -500),
          Point3D(-500, -500, -500), Point3D(500, -500, -500))
# Edges of the cube
self.a = (Segment3D(self.v[0], self.v[1]), Segment3D(self.v[1], self.v[2]),
          Segment3D(self.v[2], self.v[3]), Segment3D(self.v[3], self.v[0]),
          Segment3D(self.v[0], self.v[4]), Segment3D(self.v[1], self.v[5]),
          Segment3D(self.v[2], self.v[6]), Segment3D(self.v[3], self.v[7]),
          Segment3D(self.v[4], self.v[5]), Segment3D(self.v[5], self.v[6]),
          Segment3D(self.v[6], self.v[7]), Segment3D(self.v[7], self.v[4]))
# Algorithm for getting 4 points
def plano_limites(self, point1, point2, point3):
    plano = Plane(Point3D(point1), Point3D(point2), Point3D(point3))
    good = []
    for i in range(12):
        a = intersection(plano, self.a[i]) # Sympy intersection
        if a:
            good.append(a[0])
    return good

Answer 1

First of all, be aware that your intersection can result in more or less than four vertices. But since the region will always be convex, you can simply draw it with a triangle fan.

To sort your vertices, you need the normal n of the plane and the centroid c of all the vertices v_i:

c = 1/(number of vertices) * (v_1 + v_2 + v3 + ...)

Then, we need a coordinate system whose z-axis is the normal. For this, we can simply define an arbitrary other direction vector d and define x = normalize(cross(d, normal)), y = cross(normal, x). For cases where d conincides with normal, we need an alternative d.

We can then calculate a representative angle of any vertex in this coordinate system:

angle_i = atan2(dot(x, v_i - c), dot(y, v_i - c))

Sort by this angle and you are done.

Answer 2

Thanks to Nico's answer I was able to learn how to do this for myself but I thought it would be useful to do a full write up as it still took me further learning to understand what was going on here.

We essentiality need to find an x- and y-axis where the z-axis is aligned with the normal, we do this by taking the cross product of an arbitrary vector other. If the other vector coincides with the normal we will need to use a different vector. Two vectors coincide if their dot product is equal to 1.

By using the other vector of (0, 1, 0) (or (0, 0, -1) in the case they coincide) in a right-hand coordinate system we will produce easier to understand x- and y-axis vectors. This does not matter for the sake of the algorithm but I found this was the most crucial part initially missing in my own understanding.

By using the other vector of (0, 1, 0) when the normal of the plane is (0, 0, 1) the x-axis will be (1, 0, 0) and the y-axis will be (0, 1, 0).

By using the other vector of (0, 0, -1) when the normal of the plane is (0, 1, 0) the x-axis will be (1, 0, 0) and the y-axis will be (-1, 0, 0).

This can easily be modeled by using your right hand and turning your finger that is the normal to point up to positive z-axis (towards yourself).

def order_on_plane(vertices, normal):
    # Find the centroid of the vertices
    centroid = (0, 0, 0)
    for v in vertices:
        centroid = add(centroid, v)
    centroid = scale(centroid, 1 / len(vertices))

    # Determine the 'other' vector, used to create the axis vectors
    other = (0, 1, 0)
    
    # If the other vector coincides with the normal vector, we need to use a different other
    if math.fabs(math.fabs(dot(other, normal)) - 1.0) < 0.0001:
        other = (0, 0, -1)

    # Create the axis vectors
    x_axis = normalize(cross(other, normal))
    y_axis = normalize(cross(normal, x_axis))

    # Sort by angle as a vector from the centroid 
    angles = []
    for v in vertices:
        vector = sub(v, centroid)
        x_pos = dot(vector, x_axis)
        y_pos = dot(vector, y_axis)

        # y_pos is passed in first for east counter clockwise convention
        angle = math.atan2(y_pos, x_pos)
        angles.append(angle)

    # Sort vertices by angle
    vertices = sorted(zip(angles, vertices), key=lambda x: x[0])
    vertices = list(map(lambda x: x[1], vertices))

    return vertices