Finding rectangle vertices from a parametrized plane

I have a set of points from an .xyz file and I need to parametrize the vertices of an identified plane. I tried using the Convex Hull function but that only returned the perimeter.

My code starts like this:

#importing data
data = np.loadtxt("Asm.xyz", dtype=float, comments="L")
x = data[:, 0]
y = data[:, 1]
z = data[:, 2]

#user selects plane of interest
plane_max = float(input("Enter the largest z-value for the plane range: "))
plane_min = float(input("Enter the smallest z-value for the plane range: "))
plane = (z >= plane_min) & (z <= plane_max)

#using linear regression to parametrize plane
X = np.column_stack((x[plane], y[plane]))
Y = z[plane]

model = LinearRegression()
model.fit(X, Y)

a, b = model.coef_
c = model.intercept_

#grab points coinciding with plane
predicted_z = a * x + b * y + c
threshold = 0.001  
on_plane_indices = np.where(np.abs(predicted_z - z) < threshold)[0]
points_on_plane = data[on_plane_indices]

#attempting to use convex hull
hull = ConvexHull(points_on_plane[:,:2])
corner_points = points_on_plane[hull.vertices]

Plotting the "vertices" and my data shows that convex hull is returning the perimeter around the planes, and not just the corners like I wanted:

Convex hull and my points plotted:

I ONLY need the corners/vertices of the rectangle. My colleague suggested I use the sklearn.decomposition.PCA function but I'm not sure if that's the right way to go.

I fixed my issue by making a 2d projection of the rectangle into the XY plane:

def identify_corners(a,b,c,x,y,z,data):
    #finding which points are on the identified plane
    predicted_z = a * x + b * y + c
    threshold = 0.001
    on_plane_indices = np.where(np.abs(predicted_z - z) < threshold)[0]
    points_on_plane = data[on_plane_indices]

    #using 2d projection to find extrema of the points
    projected_points = points_on_plane[:, :2]  # do xy projection to find extrema (plane corners)
    min_x, min_y = np.min(projected_points, axis=0)
    max_x, max_y = np.max(projected_points, axis=0)
    bounding_corners_2d = np.array([
        [min_x, min_y],
        [min_x, max_y],
        [max_x, min_y],
        [max_x, max_y]])

    #finding corresponding points in case they are different in 3d space
    indices = pairwise_distances_argmin(bounding_corners_2d, projected_points) 
    corners_3d = points_on_plane[indices]
    return corners_3d