You can treat this as a clustering problem, where the cluster "centers" are actually straight lines. To compute the clustering you can use a k-means algorithm:
- Pick 4 random pairs of points. Fit a line to each, so you have 4 lines going through the point cloud.
- Repeat until the solution seems to have converged:
- For each of the points, compute the distance to each of the 4 lines.
- Assign the point to a bucket that corresponds to the line it lies closest to.
- Fit 4 new lines to the points in each of the 4 buckets (you can use linear regression or SVD)
To improve the first step you could take your idea of taking a histogram over the angles, and assign each point initially to a bucket that corresponds to the closest peak. Then fit lines to the four buckets, and start iterating.
You can also treat this as an optimization problem: pick 4 points so that the area of the difference (white area inside and black area outside of the quadrilateral) is smallest possible. Generic optimization algorithms probably work, but to make it fast you need a reasonable algorithm to compute areas.