Find the Center of the Largest Blob

Question

I have a vector of points that form blobs inside a 2d region. This 2d region represents an image. I need to find the center point of the largest blob very quickly. My algorithm that I use takes ~1 minute for a 50x50 image and it took so long that I canceled it for a 2560x1440 image. I want to be able to calculate a 2560x1440 image in just a couple seconds. Is this possible? I don't mind if the center is just an approximation. It doesn't have to be exact.

Here is an image I created from the points:

The white points are the ones that I am interested in. I need the center of the largest blob. The points that are identified could look wildly different, so take this image with a grain of salt.

I took the code from here and adapted it for my purpose. The point.h file contains an iteration value. I create a 2d array where every point outside the blobs has an iteration value of -1 and the rest have a value greater than or equal to 0 and passed that into the method.

Point.h

class Point
{
private:
    int x, y, iteration;
public:
    Point(int Xin, int Yin) : x(Xin), y(Yin), iteration(-1) {}
    Point(int Xin, int Yin, int Iterationin) : x(Xin), y(Yin), iteration(Iterationin) {}
    int getX() { return x; };
    int getY() { return y; };
    int getIteration() { return iteration; };
};

Center Point

int calculateFramePerturbation2::fill(vector<vector<int>> arr, int r, int c, vector<Point>& points) {
    int count = 0;
    if (r < arr.size() && arr[r][c] >= 0) {
        for (int i = c; i >= 0 && arr[r][i] >= 0; i--) {
            Point pt(r, i, arr[r][i]);
            points.push_back(pt);

            arr[r][i] = -1;
            count += fill(arr, r + 1, i, points) + 1;
        }
        for (int i = c + 1; i < arr[r].size() && arr[r][i] >= 0; i++) {
            Point pt(r, i, arr[r][i]);
            points.push_back(pt);

            arr[r][i] = -1;
            count += fill(arr, r + 1, i, points) + 1;
        }
    }
    return count;
}

Point calculateFramePerturbation2::print_components(vector<vector<int>> arr) {
    vector<int> pointAreas;
    vector<vector<Point>> allPoints;

    for (int r = 0; r < arr.size(); ++r) {
        for (int c = 0; c < arr[r].size(); ++c) {
            if (arr[r][c] >= 0) {
                vector<Point> points;

                fill(arr, r, c, points);
                //cout << fill(arr, r, c, points) << " ";
                //cout << points.size() << endl;
                allPoints.push_back(points);
                pointAreas.push_back(points.size());
            }
        }
    }

    //Identify the blob with the highest area.
    int max = -1;
    int index2 = 0;
    for (int i = 0; i < pointAreas.size(); i++) {
        if (pointAreas[i] > max) {
            max = pointAreas[i];
            index2 = i;
        }
    }

    cout << endl;
    cout << max << " " << index2 << endl;

    //Find the center point
    int totalX = 0;
    int totalY = 0;
    for (int i = 0; i < allPoints[index2].size(); i++) {
        totalX += allPoints[index2][i].getX();
        totalY += allPoints[index2][i].getY();
    }

    Point finalPoint(totalX / allPoints[index2].size(), totalY / allPoints[index2].size(), -1);
    return finalPoint;
}
void run(vector<Point> pointsRemaining) {
    //I'd prefer a solution where I could just take the points remaining array (it defines the blobs instead of converting it into a 2d array. 

    vector<vector<int>> glitch2DArray;

    glitch2DArray.resize(WIDTH, vector<int>(HEIGHT, -1));
    for (int i = 0; i < pointsRemaining.size(); i++) {
        glitch2DArray[pointsRemaining[i].getX()][pointsRemaining[i].getY()] = pointsRemaining[i].getIteration();
    }

    Point finalPoint = print_components(glitch2DArray);
}

Answer 1

The observation of @juvian, which pointed to the obvious mistake of copying vectors, is good. But I think that the main point is that of recognizing this is a pretty famous problem in Image Processing: Connected Components Labeling. Flood fill is one of the slowest solutions, because of its poor cache usage.

Here I'm using the SAUF algorithm published in Two Strategies to Speed up Connected Component Labeling Algorithms.

In particular this solution works by making two scans of the image and using Union-Find to manage label equivalences. One thing to note is that the second scan is not needed to solve the specific problem, because you can stop without getting the labeled image and just use the vector of blobs. This saves about 20% of the total time.

This takes few milliseconds for the specific image.

I'm sorry for posting this extremely long piece of code, but a lot of additional code is needed for having a working example. What is needed is a way to store images, the UF data structure (here everything is static, but this is not needed), the PerformLabeling function. The PNG image was converted to PBM format for easy loading.

#include <vector>
#include <fstream>
#include <string>
#include <cstdint>
#include <tuple>

using namespace std;

template <typename T>
struct mat {
    size_t rows_, cols_;
    vector<T> data_;

    mat(size_t rows = 0, size_t cols = 0) : rows_(rows), cols_(cols), data_(rows*cols) {}

    const T& operator()(size_t r, size_t c) const { return data_[cols_*r + c]; }
    T& operator()(size_t r, size_t c) { return data_[cols_*r + c]; }

    const T* operator[](size_t r) const { return &data_[cols_*r]; }
    T* operator[](size_t r) { return &data_[cols_*r]; }
};

bool loadPBM(const string& filename, mat<uint8_t>& img)
{
    ifstream is(filename, ios::binary);
    if (!is)
        return false;

    string type;
    size_t w, h;
    is >> type >> w >> h;
    if (type!="P4" || is.get()!='\n')
        return false;

    img = mat<uint8_t>(h, w);

    for (size_t r = 0; r < h; ++r) {
        uint8_t nbits = 0;
        uint8_t buffer;
        for (size_t c = 0; c < w; ++c) {
            if (nbits == 0) {
                if (!is.get(reinterpret_cast<char&>(buffer)))
                    return false;
                nbits = 8;
            }
            img(r, c) = 255 * (1 - (buffer >> --nbits & 1));
        }
    }

    return is.good();
}

bool savePGM(const string& filename, const mat<uint8_t>& img)
{
    ofstream os(filename, ios::binary);
    if (!os)
        return false;

    os << "P5\n" << img.cols_ << " " << img.rows_ << "\n255\n";
    os.write(reinterpret_cast<const char*>(&img.data_[0]), img.cols_*img.rows_);

    return os.good();
}

// Union-find (UF)
class UF {
public:
    static void Alloc(unsigned max_length) {
        P_ = new unsigned[max_length];
    }
    static void Dealloc() {
        delete[] P_;
    }
    static void Setup() {
        P_[0] = 0;// First label is for background pixels
        length_ = 1;
    }
    static unsigned NewLabel() {
        P_[length_] = length_;
        return length_++;
    }
    static unsigned GetLabel(unsigned index) {
        return P_[index];
    }
    static unsigned GetLength() {
        return length_;
    }
    static unsigned Merge(unsigned i, unsigned j)
    {
        // FindRoot(i)
        while (P_[i] < i) {
            i = P_[i];
        }
        // FindRoot(j)
        while (P_[j] < j) {
            j = P_[j];
        }
        if (i < j)
            return P_[j] = i;
        return P_[i] = j;
    }
    static unsigned Flatten()
    {
        unsigned k = 1;
        for (unsigned i = 1; i < length_; ++i) {
            if (P_[i] < i) {
                P_[i] = P_[P_[i]];
            }
            else {
                P_[i] = k;
                k = k + 1;
            }
        }
        return k;
    }
private:
    static unsigned *P_;
    static unsigned length_;
};
unsigned* UF::P_;
unsigned UF::length_;

struct blob {
    size_t x = 0, y = 0, count = 0;
    size_t left = SIZE_MAX, right = 0, top = SIZE_MAX, bottom = 0;
    blob() {}
    void accumulate(size_t r, size_t c) {
        x += c;
        y += r;
        ++count;

        if (left > c)
            left = c;
        if (right < c)
            right = c;
        if (top > r)
            top = r;
        if (bottom < r)
            bottom = r;
    }
    void accumulate(const blob& b) {
        x += b.x;
        y += b.y;
        count += b.count;

        if (left > b.left)
            left = b.left;
        if (right < b.right)
            right = b.right;
        if (top > b.top)
            top = b.top;
        if (bottom < b.bottom)
            bottom = b.bottom;
    }
};

template <typename Blob>
size_t PerformLabeling(const mat<uint8_t>& img, mat<unsigned>& img_labels, vector<Blob>& blobs)
{
    size_t h = img.rows_;
    size_t w = img.cols_;

    img_labels = mat<unsigned>(h, w); // Allocation + initialization of the output image

    UF::Alloc(((img.rows_ + 1) / 2) * ((img.cols_ + 1) / 2) + 1); // Memory allocation of the labels solver
    UF::Setup(); // Labels solver initialization

    vector<Blob> temp_blobs(((img.rows_ + 1) / 2) * ((img.cols_ + 1) / 2) + 1);

    // Rosenfeld Mask
    // +-+-+-+
    // |p|q|r|
    // +-+-+-+
    // |s|x|
    // +-+-+

    // First scan
    for (size_t r = 0; r < h; ++r) {
        // Get row pointers
        const unsigned char *img_row = img[r];
        const unsigned char *img_row_prev = img_row - img.cols_;
        unsigned *img_labels_row = img_labels[r];
        unsigned *img_labels_row_prev = img_labels_row - img_labels.cols_;

        for (size_t c = 0; c < w; ++c) {
#define CONDITION_P c > 0 && r > 0 && img_row_prev[c - 1] > 0
#define CONDITION_Q r > 0 && img_row_prev[c] > 0
#define CONDITION_R c < w - 1 && r > 0 && img_row_prev[c + 1] > 0
#define CONDITION_S c > 0 && img_row[c - 1] > 0
#define CONDITION_X img_row[c] > 0

#define ACTION_1 // nothing to do 
#define ACTION_2 img_labels_row[c] = UF::NewLabel(); // new label
#define ACTION_3 img_labels_row[c] = img_labels_row_prev[c - 1]; // x <- p
#define ACTION_4 img_labels_row[c] = img_labels_row_prev[c]; // x <- q
#define ACTION_5 img_labels_row[c] = img_labels_row_prev[c + 1]; // x <- r
#define ACTION_6 img_labels_row[c] = img_labels_row[c - 1]; // x <- s
#define ACTION_7 img_labels_row[c] = UF::Merge(img_labels_row_prev[c - 1], img_labels_row_prev[c + 1]); // x <- p + r
#define ACTION_8 img_labels_row[c] = UF::Merge(img_labels_row[c - 1], img_labels_row_prev[c + 1]); // x <- s + r

            if (CONDITION_X) {
                if (CONDITION_Q) {
                    ACTION_4
                }
                else {
                    if (CONDITION_R) {
                        if (CONDITION_P) {
                            ACTION_7
                        }
                        else {
                            if (CONDITION_S) {
                                ACTION_8
                            }
                            else {
                                ACTION_5
                            }
                        }
                    }
                    else {
                        if (CONDITION_P) {
                            ACTION_3
                        }
                        else {
                            if (CONDITION_S) {
                                ACTION_6
                            }
                            else {
                                ACTION_2
                            }
                        }
                    }
                }
            }
            else {
                ACTION_1
                continue;
            }

            temp_blobs[img_labels_row[c]].accumulate(r, c);
        }
    }

    size_t n_labels = UF::Flatten();

    blobs.resize(n_labels);
    for (size_t i = 0; i < UF::GetLength(); ++i) {
        blobs[UF::GetLabel(i)].accumulate(temp_blobs[i]);
    }

    // Second scan
    for (size_t r = 0; r < img_labels.rows_; ++r) {
        unsigned *img_row_start = img_labels[r];
        unsigned *img_row_end = img_row_start + img_labels.cols_;
        for (; img_row_start != img_row_end; ++img_row_start) {
            *img_row_start = UF::GetLabel(*img_row_start);
        }
    }

    UF::Dealloc(); // Memory deallocation of the labels solver

    return n_labels;
}

int main() {
    mat<uint8_t> img;
    if (!loadPBM("mD6sL.pbm", img))
        return EXIT_FAILURE;

    mat<unsigned> img_labels;

    vector<blob> blobs;
    auto n = PerformLabeling(img, img_labels, blobs);

    blob bmax;
    for (const auto& b :  blobs) {
        if (bmax.count < b.count)
            bmax = b;
    }
    bmax.x /= bmax.count;
    bmax.y /= bmax.count;

    mat<uint8_t> out(img.rows_, img.cols_);
    copy(begin(img_labels.data_), end(img_labels.data_), begin(out.data_));

    for (int i = -50; i <= 50; ++i) {
        out(bmax.y + i, bmax.x) = 255;
        out(bmax.y, bmax.x + i) = 255;
    }
    for (size_t i=bmax.left;i<=bmax.right;++i) {
        out(bmax.top, i) = 255;
        out(bmax.bottom, i) = 255;
    }
    for (size_t i = bmax.top; i <= bmax.bottom; ++i) {
        out(i, bmax.left) = 255;
        out(i, bmax.right) = 255;
    }

    if (!savePGM("out.pgm", out))
        return EXIT_FAILURE;

    return EXIT_SUCCESS;
}

Answer 2

Check my comment to fix your current version. Since you also asked for a version that only used your pointsRemaining, here it is. Note that this code does not compile, but should get you the idea of doing flood fill with this approach.

The worst complexity of this is O(width * height * log(width)). This might sound worse than your current version that is O(width * height), but the worst case is considering that all 2d matrix points are in pointsRemaining. A more accurate complexity of this version is O(pointsRemaining * log(width))

With this approach we avoid having to store the 2d matrix, intead just keeping the points for each row in a set for fast lookup. This approach will work very well if the amount of points from pointsRemaining is small compared to all possible points.

    void check (int y, int x, vector<set<int>> &rows, queue<pair<int, int>> &q ) {
        if (rows[y].includes(x)) {
            q.push({y, x});
            rows[y].remove(x);
        }    
    }

    void run(vector<Point> pointsRemaining) {
        vector<set<int>> rows(HEIGHT);

        for (int i = 0; i < HEIGHT; i++) {
            rows[i] = set<int>();        
        }

        for (Point p: pointsRemaining) {
            rows[p.getY()].add(p.getX());
        }

        queue<pair<int, int>> q;

        for (Point p: pointsRemaining) {
            if (rows[p.getY()].includes(p.getX())) {
                q.push({p.getY(), p.getX()});
                rows[p.getY()].remove(p.getX());
                int size = 0; 

                while (!q.empty()) { // queue based flood fill to avoid recursion depth limits
                    y, x = q.front();
                    q.pop();
                    size++;

                    check(y, x + 1, rows, q);
                    check(y, x - 1, rows, q);
                    check(y + 1, x, rows, q);
                    check(y - 1, x, rows, q);

                }

                cout << size << "\n"; // size of component (blob)
            }
        }

    }