Home > Back-end >  What is the fastest way to find all squares in an array of points?
What is the fastest way to find all squares in an array of points?

Time:12-16

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

CodePudding user response:

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

enter image description here

CodePudding user response:

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 n3.

CodePudding user response:

Similar to Ted's answer:

  • enter all points to a hash set;

  • double loop on all "first quadrant" point pairs (X1 ≥ X0 and Y1 > Y0);

    • test presence of (X0 Y0 - Y1, Y0 X1 - X0) and (X1 Y0 - Y1, Y1 X1 - X0).

The "first quadrant" condition ensures that no square will be duplicated (a square has exactly one side in each quadrant).

This takes time O(n²), which is worst-case optimal, as there can be O(n²) squares.

  • Related