既然你提到它是一个 “计算机图形坐标系（没有负坐标）"，我也假设整数值，这使得它相当直接。

• 创建Point的std::unordered_set
• 对于两个点A和B的每个组合，计算第三个点C的位置，第三个点C相对于前两个点A -> B旋转90°。
• 如果该集合中存在这样的第三个点，则计算最后一个点D相对于B -> C的位置。
• 如果找到匹配，则从集合中删除这四个点。
• 继续，直到你通过所有的点。

示例：
首先是Point类模板：

template <class T = unsigned long long> // or some other integer type
struct Point {
    bool operator==(const Point& o) const { return x == o.x && y == o.y; }
    bool operator!=(const Point& o) const { return !(*this == o); }

    Point& operator-=(const Point& o) {
        x -= o.x;
        y -= o.y;
        return *this;
    }

    T x, y;
};

template <class T>
Point<T> operator-(const Point<T>& lhs, const Point<T>& rhs) {
    Point rv = lhs;
    rv -= rhs;
    return rv;
}

template <class T>
std::ostream& operator<<(std::ostream& os, const Point<T>& p) {
    return os << '{' << p.x << ',' << p.y << '}';
}

template <class T>
struct std::hash<Point<T>> {
    std::size_t operator()(const Point<T>& p) const {
        std::hash<T> h;
        // boost::hash_combine:
        std::size_t rv = h(p.x);
        rv ^= h(p.y) + 0x9e3779b9 + (rv << 6) + (rv >> 2);
        return rv;
    }
};

这样我们就可以通过B - A得到x和y的不同之处，并且专用的std::hash可以让我们把它存储在像std::unordered_set这样的容器中，std::unordered_set要求 Key 是可哈希的。
接下来，使用一个函数模板，通过将Point<>相对于其他两个Point<>旋转90°来获得Point<>：

template<class T>
Point<T> rot90(const Point<T>& A, const Point<T>& B) {
    // C.x = B.x + (B-A).y
    // C.y = B.y - (B-A).x
    auto BA = B - A;
    return {B.x + BA.y, B.y - BA.x};
}

然后，实际的匹配过程可能如下所示：

int main() {
    std::unordered_set<Point<>> points{ /* ... */ };

    // first, second, third and last iterator:
    // fit(A), sit(B), tit(C), lit(D)
    for (auto fit = points.begin(); fit != points.end();) {
        bool found = false;

        // temporarily extract the node from the set to loop over the rest
        // of the nodes:
        auto next = std::next(fit);
        auto node = points.extract(fit);
        
        // loop over the rest of the nodes:
        auto sit = points.begin();
        decltype(sit) tit, lit;

        for(; sit != points.end(); ++sit) {
            // calculate where C should be:
            auto candidate = rot90(node.value(), *sit);
            if(tit = points.find(candidate); tit != points.end()) {
                // third Point found - try for the last too:
                candidate = rot90(*sit, candidate);
                if(lit = points.find(candidate); lit != points.end()) {
                    found = true;
                    break;
                }
            }
        }
        if(found) {
            std::cout << "FOUND: " << node.value() << ' ' << *sit << ' '
                                   << *tit << ' ' << *lit << '\n';
            // erase the points from the set
            if(next == sit) ++next; // next being erased, step next
            points.erase(sit);
            if(next == tit) ++next; // next being erased, step next
            points.erase(tit);
            if(next == lit) ++next; // next being erased, step next
            points.erase(lit);
        } else {
            // reinsert the first Point node since no square was found
            points.insert(fit, std::move(node));
        }
        fit = next; // try next
    }
}

Demo

注：

对于某些点的组合，您可能会找到比预期更少的方块。例如，上述算法（以及您的原始算法）可能会在此处找到1 * 或 * 2个方块：

{100,1} {101,1} {101,0} {100,0}
{102,1} {103,1} {103,0} {102,0}

那是因为

{101,0} {101,1} {102,1} {102,0}

也形成了一个（第三个）正方形，由于您删除了已经使用的点，如果上面的点先找到，其他两个点将不会显示。如果您希望显示所有三个点，您可以将它们收集到另一个集合中。首先是Square类模板：

template<class T = unsigned long long>
class Square {
public:
    // this requires the points to be added in the order you'd "paint" them,
    // that is, in the order they were found by the matching algorithm below
    Square(Point<T> A, Point<T> B, Point<T> C, Point<T> D) : points{A, B, C, D} {
        // rotate to put the "smallest" Point first
        auto mit = std::min_element(points.begin(), points.end());
        std::rotate(points.begin(), mit, points.end());
    }

    bool operator<(const Square& o) const {
        return points < o.points;
    }

    friend std::ostream& operator<<(std::ostream& os, const Square<T> &s) {
        return os << '{' << s.points[0] << ',' << s.points[1] << ',' 
                         << s.points[2] << ',' << s.points[3] << '}';
    }

private:
    std::array<Point<T>, 4> points;
};

匹配变得更加直接：

template<class T>
auto FindSquares(const std::unordered_set<Point<T>>& points,
                 std::set<Square<T>>& sq_found)
{
    //auto start = std::chrono::steady_clock::now();

    // loop over all points
    for (auto fit = points.begin(), end = points.end(); fit != end; ++fit) {
        // loop over all the other points
        for(auto sit = points.begin(); sit != end; ++sit) {
            if(sit == fit) continue;

            // try to find the third point:
            if(auto tit = points.find(rot90(*fit, *sit)); tit != points.end())
            {
                // third Point found - try for the last too:
                if(auto lit = points.find(rot90(*sit, *tit)); lit != points.end())
                {
                    // store this square (if it's not there already)
                    sq_found.emplace(*fit, *sit, *tit, *lit);
                }
            }
        }
    }   
    //return std::chrono::steady_clock::now() - start;
}

int main() {
    std::unordered_set<Point<>> points{/* ... */};
    std::set<Square<>> sq_found;

    FindSquares(points, sq_found);

    for(auto& sq : sq_found) {
        std::cout << "FOUND: " << sq << '\n';
    }
}

Demo
我对上面的FindSquares函数做了一些测量，在随机位置放置点，持续时间只是为了显示它的缩放方式：
| 尺寸|点数|找到的方块|持续时间|
| - ------|- ------|- ------|- ------|
| 590x591|一千五百|1-6| 0.04s |
| 590x591|一万五千|小行星35| 5.65s |
| 590x591|三万|小行星565800| 31.00s |
| 小行星1777 x1057|一千五百|无| 0.04s |
| 小行星1777 x1057|一万五千|一千| 5.40s |
| 小行星1777 x1057|三万|小行星157| 27.70s |
x1c 0d1x的数据

你必须测试两个初始点的任意组合。我们称之为A和B。但是由于你只需要每个可能的正方形一次，你可以要求第三个点，我们称之为C，以验证（AC）在三角学意义上是（AB）旋转90°。如果一个可能的候选者可以通过旋转相反的方向获得，你将从不同的初始对中找到最终的正方形。
但这里的条件非常简单：

• yC - yA = xB - xA且
• xC -xA = -（yB - yA）

当你有一个可能的候选点C，点D将必须验证（CD）=（AB）的意思是：

• xD - xC = xB - xA以及
• yD - yC = yB - yA

但是如果你使用浮点值作为坐标，你可能会遇到精度问题。我们都知道浮点运算是不精确的。所以你应该定义一个 epsilon 值（通常10-6是一个可接受的值），并用接近以下代码替换所有的相等性测试：

if (abs(x - y) < epsilon) ...

我目前无法编写和测试代码（无法访问可接受的编译器），但此算法应该可以稍微缩短时间：操作简单得多，并且由于如果你已经找到可接受的第三个点，则你仅搜索第四个点，因此复杂度应当接近n3。