Which CV algorithm should I choose to detect the corners? (Corner detection, OpenCV)

I'm working on a CV tool for image analysis and use 4 points to form a transformation matrix (map) for image correction (perspective projection - cv.four_point_transform). I have a set of source images obtained from a thermal imager with the presence of distortion. I'm successfully applying a pipeline of standard OpenCV functions. You can see it in the pictures below (Fig. 1). Please note that the main pipeline consists of the following steps:

1. Correction of distortion.
2. Smoothing with a bilateral filter.
3. Threshold.
4. Harris corner detector.
5. Erosion.
6. The weighted average value of a point on a plane over a point cloud after erosion. In the "Target" image you see 4 points, the actual size of 1 pixel (enlarged for clarity).

• Which approach do you think is easier to implement in the future?
• How to properly approach the removal of parasitic points?
• Are there simpler approaches to detection corners?

Unfortunately, I come across cases where the Harris corner detector doesn't cope and doesn't detect obtuse angles. I started testing different approaches such as:

1. Threshold -> Contours -> Approximate Contour -> Points.
2. Threshold -> Canny Edges Detection -> Dilation -> FAST.
3. Threshold -> Canny Edges Detection -> Dilation -> SIFT.
4. Threshold -> Canny Edges Detection -> Dilation -> Probabilistic Hough lines -> Bentley-Ottmann points.

As you can see, some of the approaches can be used, but they form a number of parasitic points. I'm afraid that it will be more difficult to fight them than it seems at first glance (Fig. 2).

Fig. 1 - Successful detection

Fig. 2 - Unsuccessful detection

Not so long ago, Dev Aggarwal (@dev-aggarwal) published an interesting and simple approach for determining a fitted quadrilateral. This approach works for me. See the original post at the link.

The theoretical layout of the algorithm is presented in the article - View Frustum Optimization To Maximize Object’s Image Area.

I duplicate the source code here:

import sympy

def appx_best_fit_ngon(img, n: int = 4) -> list[(int, int)]:
    # convex hull of the input mask
    img_gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
    contours, _ = cv2.findContours(
        img_gray, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE
    )
    hull = cv2.convexHull(contours[0])
    hull = np.array(hull).reshape((len(hull), 2))

    # to sympy land
    hull = [sympy.Point(*pt) for pt in hull]

    # run until we cut down to n vertices
    while len(hull) > n:
        best_candidate = None

        # for all edges in hull ( <edge_idx_1>, <edge_idx_2> ) ->
        for edge_idx_1 in range(len(hull)):
            edge_idx_2 = (edge_idx_1 + 1) % len(hull)

            adj_idx_1 = (edge_idx_1 - 1) % len(hull)
            adj_idx_2 = (edge_idx_1 + 2) % len(hull)

            edge_pt_1 = sympy.Point(*hull[edge_idx_1])
            edge_pt_2 = sympy.Point(*hull[edge_idx_2])
            adj_pt_1 = sympy.Point(*hull[adj_idx_1])
            adj_pt_2 = sympy.Point(*hull[adj_idx_2])

            subpoly = sympy.Polygon(adj_pt_1, edge_pt_1, edge_pt_2, adj_pt_2)
            angle1 = subpoly.angles[edge_pt_1]
            angle2 = subpoly.angles[edge_pt_2]

            # we need to first make sure that the sum of the interior angles the edge
            # makes with the two adjacent edges is more than 180°
            if sympy.N(angle1 + angle2) <= sympy.pi:
                continue

            # find the new vertex if we delete this edge
            adj_edge_1 = sympy.Line(adj_pt_1, edge_pt_1)
            adj_edge_2 = sympy.Line(edge_pt_2, adj_pt_2)
            intersect = adj_edge_1.intersection(adj_edge_2)[0]

            # the area of the triangle we'll be adding
            area = sympy.N(sympy.Triangle(edge_pt_1, intersect, edge_pt_2).area)
            # should be the lowest
            if best_candidate and best_candidate[1] < area:
                continue

            # delete the edge and add the intersection of adjacent edges to the hull
            better_hull = list(hull)
            better_hull[edge_idx_1] = intersect
            del better_hull[edge_idx_2]
            best_candidate = (better_hull, area)

        if not best_candidate:
            raise ValueError("Could not find the best fit n-gon!")

        hull = best_candidate[0]

    # back to python land
    hull = [(int(x), int(y)) for x, y in hull]

    return hull

Use as follows:

hull = appx_best_fit_ngon(img)

# add lines
for idx in range(len(hull)):
    next_idx = (idx + 1) % len(hull)
    cv2.line(img, hull[idx], hull[next_idx], (0, 255, 0), 1)

# add point markers
for pt in hull:
    cv2.circle(img, pt, 2, (255, 0, 0), 2)

Of the minuses, I see only one additional dependency on the SymPy package, but this is not critical.