Centroid of only specific (dark) pixels in an 8-bit image

[VeMo]

Hello,

I have a small program that calculates centroid as attached. But this program currently calculates centroid of whole image and I want to calculate centroid of only dark pixels and ignore everything else. For example, the attached image has a cross mark made with dark gray/black pixels in form of circles. I want the program to find centroid of only dark pixels and ignore whole white background. I need some kind of thresholding/restriction that will allow only pixels with certain value (for e.g. <100 because it is an 8-bit image with values of pixels between 0 and 255) or less pass through the vi that calculates centroid. Can someone please help as I am new with LabVIEW and I am trying to learn?

Thanks

VeMo 

[Andrey_Dmitriev]

There are different ways how to get this.

An easiest way probably is just invert image. The centroid computed based on intensities, so on the inverted he will sit much closer to the "true" center.

Another way to use threshold, in your case automated threshold will work quite well, so you can perform centroid computation on the thresholded image:

Another way is to analyse bound rects of the particles after thresholding and do it without centroid at all:

Then you can choose which way is most suitable for you:

Attached example, downgraded to LV 2018.

[VeMo]

Hi Andrey,

Thank you so much for your reply. I was working on making a simple manual program on LabVIEW to calculate centroid rather than using a given vi in the system so I can learn how the centroid is calculated and if it is right. Attached is the program I made for trial and I was wondering if any one of your suggestions will work for my manual centroid calculation program? I will try doing some more work but since I am new I may be little slow in updating my program. 

[Andrey_Dmitriev]

---

@VeMo wrote:

I was wondering if any one of your suggestions will work for my manual centroid calculation program? I will try doing some more work but since I am new I may be little slow in updating my program. 

---

Well, it is a classical centroid computation. Your version is OK and seems to be correct. It is not the fastest in terms of performance (it can be improved by a factor of 10x, I guess), but it's acceptable. My recommendation is to always stay away from converting IMAQ Images to Arrays if you can, because there are a lot of penalties and memory waste. Additionally, operations on native LabVIEW arrays are not very fast compared to highly optimized image processing libraries (usually vectorized). Is there a requirement not to use IMAQ Centroid, or was it just a learning exercise?

[VeMo]

Yes I am working on something and I need to have information on what the process is doing so using a direct vi is not the option right now. I am using the centroid vi just to verify my centroid calculation but it is preferred to make my own program to calculate centroid. I would appreciate if you have any other better suggestion than using an array for manual centroid calculations as you mentioned? 

[Andrey_Dmitriev]

---

@VeMo wrote:

Yes I am working on something and I need to have information on what the process is doing so using a direct vi is not the option right now. I am using the centroid vi just to verify my centroid calculation but it is preferred to make my own program to calculate centroid. I would appreciate if you have any other better suggestion than using an array for manual centroid calculations as you mentioned? 

---

OK, the first suggestion is about arrays. Even if you working with arrays, then the very first thing which I would recommend is to move from floating-point to integer arithmetic. For centroid calculation, you only need to add and multiply in the loops. But take note that you can easily overflow your accumulators, so here you will need the U64 type. This is another argument to have integers instead of floats, is because at the end you will deal with small numbers (pixel intensities) with very large numbers in shift registers and this is not good (from precision point of view) for float point arithmetic.

I will add a randomly generated image with some gradient as well as a simple time check, this is how it looks:

The next suggestion below may be slightly ahead of your knowledge, but for educational purposes, I'll do this, could be useful for some else as well.

If you take a look inside if the IMAQ Centroid, then you will see a DLL call behind the scenes:

Usually, a well-written DLL will in most cases beat native LabVIEW code. Let me illustrate this.

To get this I will need to have a pointer to the first pixel and geometry (refer to IMAQ GetImagePixelPtr VI.)

But trivial DLL will not give much, so I will use AVX2 extension, which is available on most of CPU nowadays, as result computation on multiple pixels (8) with single command will be possible and this will increase performance, so my nested for loops will be like this:

An easiest way to program this is using Intel Intrinsics, there are only fifty lines of the C code:

Spoiler

#include "framework.h"
#include <immintrin.h>
#include <inttypes.h>
#include "CentroidAVX2DLL.h"

typedef struct {
    float x;
    float y;
} Centroid;

CENTROIDAVX2DLL_API int centroid_avx2(const uint8_t* image, int width, int height, Centroid* centroid)
{
    __m256i vx_base = _mm256_set_epi32(7, 6, 5, 4, 3, 2, 1, 0);
    __m256i vsum = _mm256_setzero_si256();
    __m256i vsum_x = _mm256_setzero_si256();
    __m256i vsum_y = _mm256_setzero_si256();

    for (int y = 0; y < height; y++) {
        __m256i vx = vx_base;
        for (int x = 0; x < width; x +=  {
            __m256i vpixels = _mm256_cvtepu8_epi32(_mm_loadl_epi64((__m128i*)(image + y * width + x)));
            vsum = _mm256_add_epi64(vsum, _mm256_cvtepi32_epi64(_mm256_extracti128_si256(vpixels, 0)));
            vsum = _mm256_add_epi64(vsum, _mm256_cvtepi32_epi64(_mm256_extracti128_si256(vpixels, 1)));

            __m256i vx_pixels = _mm256_mullo_epi32(vpixels, vx);
            vsum_x = _mm256_add_epi64(vsum_x, _mm256_cvtepi32_epi64(_mm256_extracti128_si256(vx_pixels, 0)));
            vsum_x = _mm256_add_epi64(vsum_x, _mm256_cvtepi32_epi64(_mm256_extracti128_si256(vx_pixels, 1)));

            __m256i vy_pixels = _mm256_mullo_epi32(vpixels, _mm256_set1_epi32(y));
            vsum_y = _mm256_add_epi64(vsum_y, _mm256_cvtepi32_epi64(_mm256_extracti128_si256(vy_pixels, 0)));
            vsum_y = _mm256_add_epi64(vsum_y, _mm256_cvtepi32_epi64(_mm256_extracti128_si256(vy_pixels, 1)));

            vx = _mm256_add_epi32(vx, _mm256_set1_epi32(8));
        }
    }

    uint64_t sum[4], sum_x[4], sum_y[4];
    _mm256_storeu_si256((__m256i*)sum, vsum);
    _mm256_storeu_si256((__m256i*)sum_x, vsum_x);
    _mm256_storeu_si256((__m256i*)sum_y, vsum_y);

    uint64_t total_sum = sum[0] + sum[1] + sum[2] + sum[3];
    uint64_t total_sum_x = sum_x[0] + sum_x[1] + sum_x[2] + sum_x[3];
    uint64_t total_sum_y = sum_y[0] + sum_y[1] + sum_y[2] + sum_y[3];

    centroid->x = (float)((double)total_sum_x / total_sum);
    centroid->y = (float)((double)total_sum_y / total_sum);

    return 0;
}

Here how this DLL will be called (it almost the same like NI does in IMAQ GetImagePixelPtr Example):

Now its a time to put all together and perform simple benchmark and test (the results must be obviously the same from all three algorithms), I'll try to make as compact as possible, you can do it with SubVIs, of course:

Above all three algorithm called sequentially, the time of the execution is measured. The size of the image was increased to 8K to get meaningful time.

Results:

As you can see, AVX2 - based DLL is around 20x faster than LabVIEW and 4x faster than IMAQ Centroid, and the results are exactly the same:

Full source code in the attachment, downgraded to LV2019; Visual Studio 2022 v.17.12.1 was used for DLL. This code is x64 only.