What is the fastest way to find all squares in an array of points?

Question

The main problem of my program is speed. I mean it works, but terribly slowly.

I have an array of points and I need to find squares here, so in order to do this I need to check all possible combinations of 4 points. So I used three nested loops. In general it would take about n⁴ operations.

If I had up to 100 elements in the array it would be more or less normal, but I have 500 elements, so 500⁴, which is quite a lot. Is here any faster way to do this?

/*j, o, k, l are numbers of each point. This loop will check all possible combinations
    (except 2 or more points have the  same number(it means also the same coordinades) or they already are part of some figure*/ 
    for (int j=0; j < i-1; j++)
    {      
        if (Point[j].p != true)
        {
            for (int o = 1; o < i - 2; o++)
            {
                if ((o != j) && (Point[o].p != true))
                {
                    for (int k = 2; k < i - 3; k++)
                    {
                        if ((k!= j) && (k != o) && (Point[k].p != true))
                        {
                            for (int l = 3; l < i - 4; l++)
                            {
                                if ((l != k) && (Point[l].p != true) && (l != o) && (l!=j))
                                {
                                    vx1 = abs(Point[o].x - Point[j].x); //vectors coordinates
                                    vx2 = abs(Point[k].x - Point[o].x);
                                    vy1 = abs(Point[o].y - Point[j].y);
                                    vy2 = abs(Point[k].y - Point[o].y);
                                    vx3 = abs(Point[l].x - Point[k].x);
                                    vy3 = abs(Point[l].y - Point[k].y);
                                    vx4 = abs(Point[j].x - Point[l].x);
                                    vy4 = abs(Point[j].y - Point[l].y);
                                    dx1 = abs(Point[k].x - Point[j].x); //diagonals coordinates
                                    dy1 = abs(Point[k].y - Point[j].y);
                                    dx2 = abs(Point[o].x - Point[l].x);
                                    dy2 = abs(Point[o].y - Point[l].y);
                                    v1 = sqrt(vx1 * vx1 + vy1 * vy1); // sides length
                                    v2 = sqrt(vx2 * vx2 + vy2 * vy2);
                                    v3 = sqrt(vx3 * vx3 + vy3 * vy3);
                                    v4 = sqrt(vx4 * vx4 + vy4 * vy4);
                                    d1 = sqrt(dx1 *dx1 + dy1 * dy1); //diagonals length
                                    d2 = sqrt(dx2 * dx2 + dy2 * dy2);
                                    if (
                                        ((abs(v1-v2)<=0.5) && (abs(v3-v4)<=0.5) && (abs(v3-v2)<=0.5) && (v1<d1)) /*cheks all  sides are equal and if the diagonal is bigger than side*/  
                                        && (Point[k].p != true && Point[o].p != true && Point[j].p != true)/*checks if the points aren`t the part of any figure yet*/ 
                                        &&(abs(d1 - d2)<=0.5)/*checks if the diagonals are equal*/)
                                    {
                                        q++;
                                        Point[j].p = true; // it means that point with this number is already a part of some figure, so it will be skipped in next operations
                                        Point[o].p = true;
                                        Point[k].p = true;
                                        Point[l].p = true;
                                        // the output
                                        out << "Figure " << q << ":" << "x1=" << Point[j].x << " y1=" << Point[j].y << " x2="
                                            << Point[o].x << " y2=" << Point[o].y <<
                                            " x3=" << Point[k].x << " y3=" << Point[k].y <<
                                            " x4=" << Point[l].x << " y4=" << Point[l].y << endl;
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }

Answer 1

Since you mentioned it's a "computer graphics coordinate system(without negative coordinates)" I'm also assuming integer values which makes it rather straight forward.

• Create a std::unordered_set of Point
• For each combination of two points, A and B, calculate the position of a third point, C, rotated 90° in relation to the first two, A -> B.
• If such a third point exists in the set, calculate the position of the last point, D, in relation to B -> C.
• If a match is found, remove the four points from the set.
• Continue until you've gone through all points.

Example:

First, a Point class template:

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;
    }
};

This will let us do B - A to get a Point with the x and y difference and the specialized std::hash will let us store it in containers like std::unordered_set that requires the Key to be hashable. The operator<< overload is just to make printing Point<> instances easier.

Next, a function template to get a Point<> from rotating 90° in relation to two other Point<>s:

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};
}

Then doing the actual matching could look like this:

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

---

Note:

For certain combinations of points you may find fewer squares than expected. For example, the above algorithm (as well as your original algorithm) may find 1 or 2 squares here:

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

That's because the points

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

also forms a (third) square and since you remove the points already used, the other two will not show up if the above is found first. If you instead want all to show all three you can collect them in another set. I'll use a std::set to store the squares which requires operator< for the Square objects. To simplify the implementation of the Square::operator<, first add Point::operator<:

    bool operator<(const Point& o) const {
        return std::tie(x, y) < std::tie(o.x, o.y);
    }

Then, the Square class template:

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; // uses Point::operator< 
    }

    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;
};

And the matching becomes somewhat even more straight forward:

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

---

I made a few measurements of the FindSquares function above when placing points at random places. The durations are only meant to show how it scales:

Dimensions	Points	Squares found	Duration
590x591	1500	1-6	0.04s
590x591	15000	35050	5.65s
590x591	30000	565800	31.00s
1777x1057	1500	0	0.04s
1777x1057	15000	1000	5.40s
1777x1057	30000	15700	27.70s

Answer 2

You will have to test for any combination of 2 initial points. Let us call them A and B. But as you only want each possible square once, you can require the third point, let us call it C, to verify (AC) is (AB) turned with 90° in trigonometric sense. If a possible candidate could be obtained by turning in the opposite sense, you would find the final square from a different initial pair.

But here the condition is damned simple:

• yC - yA = xB - xA and
• xC -xA = - (yB - yA)

And when you have a possible candidate for point C, point D will have to verify (CD) = (AB) meaning:

• xD - xC = xB - xA and
• yD - yC = yB - yA

But if you use floating point values for your coordinates, you could fall in an accuracy problem. We all know that floating point arithmetics is broken inaccurate. So you should define an epsilon value (generally 10^-6 is an acceptable value) and replace all equality tests with code close to:

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

I cannot currently write and test code (no access to a acceptable compiler) but this algorithm should slightly reduce time: operations are much simpler and as you only search the 4th point if you have found an acceptable 3rd one the complexity should be closer to n³.