Implementing LSL for Connected Component Labeling/Blob Extraction

Question

I am trying to implement the Light Speed Labeling algorithm from https://www.lri.fr/~lacas/Publications/JRTIP10.pdf. I tried to follow the algorithmic approach described in the paper as closely as possible (pgs 9-12) but the output doesn't make sense after the equivalence stage.

Anyone have any idea what the problem is?

void segment(const unsigned * Xi, const unsigned& N, unsigned * ERi, unsigned *RLCi, unsigned& ner) 
{ 
unsigned x1     = 0;
unsigned f      = 0;
unsigned er     = 0;    

for (unsigned j = 0; j < N; ++j)
{
    const unsigned x0 = Xi[j];
    f = x0 ^ x1;
    RLCi[er] = j;         
    er = er + f;
    ERi[j] = er;         
    x1 = x0;        
}    
if (x1 != 0)
{
    RLCi[er] = N;
}
er = er + x1;
ner = er;
}

void equivalance(const unsigned& ner, unsigned * RLCi, unsigned * EQ, unsigned * ER0, unsigned * ERA0, unsigned * ERA1, unsigned& nea, const unsigned& N = 0)
{
for (unsigned er = 1; er <= ner; er += 2)
{        
    int j0 = RLCi[er - 1];
    int j1 = RLCi[er] - 1;
   // Unnecessary given optimization and need for 4-connectivity:
    // if (j0 > 0) j0 = j0 - 1;
   // if (j1 < N - 1) j1 = j1 + 1;

    int er0 = ER0[j0];
    int er1 = ER0[j1];

    if (!(er0 & 1)) er0 = er0 + 1;
    if (!(er1 & 1)) er1 = er1 - 1;       

    if (er1 >= er0) // adjacent label
    {            
        unsigned ea = ERA0[er0];
        unsigned a = EQ[ea];
        for (unsigned erk = er0 + 2; erk <= er1; ++erk)
        {
            unsigned eak = ERA0[erk];              
            unsigned ak = EQ[eak];
            if (a < ak)
            {
                EQ[eak] = a;
            }
            else
            {
                a = ak;                    
                EQ[ea] = a;
                ea = eak;
            }                
        }           
        ERA1[er] = a;
    }
    else
    {
        nea = nea + 1;
        ERA1[er] = nea;
    }
}   
}

typedef std::vector<unsigned> value_type;
void bwlabel(const double* X, unsigned * EA, const unsigned& N, const unsigned& M)
{
unsigned nea = 0;
const unsigned size = N * M;

value_type EQ(size, 0), ER(size, 0), ERA(size, 0), A(size, 0), RLC(M * (2 * N), 0), IN(X, X + size), NER(M, 0);

// Step 1
for (int m = 0; m < M; ++m)
{        
    segment(&IN[0] + N * m, N, &ER[0] + N * m, &RLC[0] + m * (2 * N), NER[m]);       
}    
// Step 2
for (int m = 1; m < M; ++m)
{        
    equivalance(NER[m], &RLC[0] + m * (2 * N), &EQ[0], &ER[0] + (m - 1) * N, &ERA[0] + (m - 1) * N, &ERA[0] + m * N, nea, N);               
}
// Step 3
for (int j = 0; j < size; ++j)
{
    EA[j] = ERA[ER[j]];        
} 

// Step 4
unsigned na = 0;
for (int e = 0; e < size; ++e)
{
    if (EQ[e] != e)
    {
        A[e] = EQ[EQ[e]];
    }
    else
    {
        na = na + 1;
        A[e] = na;
    }
}
// Step 5
for (int j = 0; j < size; ++j)
{
    EA[j] = A[EA[j]];
}
}
IN=
1 1 0 1 1 0 0 1 1 1 
0 1 1 0 1 0 0 1 1 1 
1 0 1 1 1 1 1 0 1 0 
1 0 0 0 0 1 1 0 1 0 
0 0 1 1 0 0 0 1 1 1 
0 1 1 0 0 0 1 0 1 0 
1 0 1 1 1 1 1 0 0 0 
1 0 1 0 1 0 0 0 1 0 
0 1 1 1 1 0 1 1 0 1 
0 0 1 1 1 0 1 0 0 0 


RLC=
[0]0 [1]2 [2]3 [3]5 [4]7 [5]10 [6]0 [7]0 [8]0 [9]0 
[0]0 [1]0 [2]0 [3]0 [4]0 [5]0 [6]0 [7]0 [8]0 [9]0 
[0]1 [1]3 [2]4 [3]5 [4]7 [5]10 [6]0 [7]0 [8]0 [9]0 
[0]0 [1]0 [2]0 [3]0 [4]0 [5]0 [6]0 [7]0 [8]0 [9]0 
[0]0 [1]1 [2]2 [3]7 [4]8 [5]9 [6]0 [7]0 [8]0 [9]0 
[0]0 [1]0 [2]0 [3]0 [4]0 [5]0 [6]0 [7]0 [8]0 [9]0 
[0]0 [1]1 [2]5 [3]7 [4]8 [5]9 [6]0 [7]0 [8]0 [9]0 
[0]0 [1]0 [2]0 [3]0 [4]0 [5]0 [6]0 [7]0 [8]0 [9]0 
[0]2 [1]4 [2]7 [3]10 [4]0 [5]0 [6]0 [7]0 [8]0 [9]0 
[0]0 [1]0 [2]0 [3]0 [4]0 [5]0 [6]0 [7]0 [8]0 [9]0 


ER=
1 1 2 3 3 4 4 5 5 5 
0 1 1 2 3 4 4 5 5 5 
1 2 3 3 3 3 3 4 5 6 
1 2 2 2 2 3 3 4 5 6 
0 0 1 1 2 2 2 3 3 3 
0 1 1 2 2 2 3 4 5 6 
1 2 3 3 3 3 3 4 4 4 
1 2 3 4 5 6 6 6 7 8 
0 1 1 1 1 2 3 3 4 5 
0 0 1 1 1 2 3 4 4 4 


ERA=
0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 
0 1 0 0 0 0 0 0 0 0 
0 0 0 2 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 3 0 0 
0 0 0 4 0 5 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 


EQ=
0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 


EA=
1 1 1 1 1 1 1 1 1 1 
1 1 1 1 1 1 1 1 1 1 
1 1 1 1 1 1 1 1 1 1 
1 1 1 1 1 1 1 1 1 1 
1 1 1 1 1 1 1 1 1 1 
1 1 1 1 1 1 1 1 1 1 
1 1 1 1 1 1 1 1 1 1 
1 1 1 1 1 1 1 1 1 1 
1 1 1 1 1 1 1 1 1 1 
1 1 1 1 1 1 1 1 1 1 

Answer 1

For starters, you have omitted the right border compensation variable b. This affects the results of RLC.

You are obtaining (for the first row):

[0]0 [1]2 [2]3 [3]5 [4]7 [5]10 [6]0 [7]0 [8]0 [9]0

when you should be obtaining:

[0]0 [1]1 [2]3 [3]4 [4]7 [5]10 [6]0 [7]0 [8]0 [9]0

Thanks.

Answer 2

I have tried to follow the same article. The article has multiple errors and misleading parts.

EQ for example, is a Union-Find data structure, not a plain array.

Here is a working C implementation I wrote based on YACCLAB version of LSL (https://github.com/prittt/YACCLAB/blob/master/include/labeling_lacassagne_2016_code.inc):

// input type
typedef int Xt;
// type for tables, change to unsigned short 
// if you won't process images larger than 256x256 
typedef unsigned uint;

// tables used in LSL:
// X  [H][W];               // input b&w image
// A  [H][W];               // output image (labels)
// ER [H][W];
// ERA[H][W + 1];
// RLC[H][(W + 1) & ~1];    // (RLC width is rounded to next multiple of 2)
// NER[H];

// YACCLAB uses "compact RLC" - an array of pointers to RLC table's rows.
// my benchmark shows that using it does not improve performance,
// while consuming additional memory for it's pointer array.
// uncomment lines related to 'NRLC' to enable.
// NRLC[H];

// specify pointer to each table manually,
// (you might want to keep RLC and deallocate everything else)
// or allocate everything in a single block and use 'lsl_tables_from_ptr()'.
// no initialization of allocated memory needed.
struct lsl_t {
    uint *ER, *ERA, *RLC, *NER, *EQ; //, *NRLC;
    size_t W, H;
};

// total size of all tables (in integers)
size_t lsl_tables_size(size_t W, size_t H);

//
lsl_t lsl_tables_from_ptr(void* memory, size_t W, size_t H);

// lsl
size_t lsl(Xt const in[], uint out[], lsl_t const* a);

// a union-find data structure
struct uf_t {
    uint* EQ;   // label equivalence table
    size_t n;   // number of equivalence records
};

//
uint eq_new_label(uf_t* uf);
// get label that 'i' is equivalent to
uint eq_get_label(uf_t const* uf, uint i);
// 'root' is a "true" or a "final" label,
// a label that is not equivalent to any other labels
uint eq_find_root(uf_t const* uf, uint er);
// set 'e' as equivalent to 'r'
void eq_update_table(uf_t* uf, uint e, uint r);
// collapses long equivalence chains
uint eq_flatten(uf_t* uf);

// this algorithm detects "segments"
// in each independent row of the input image.
// a "segment" is a span of consecutive
// background or foreground pixels.
// background pixels are labeled with even numbers
// and foreground pixels with odd.
void seg(
    Xt const Xi[], uint w,
    uint ERi[], uint RLCi[], uint* ner)
{
    uint er = 0;
    int x0, x1 = 0; // curr and prev val of X
    int b = 0;      // right border compensation
    int f = 0;      // front detection

    for (size_t j = 0; j < w; ++j) {
        x0 = Xi[j];
        f = x0 ^ x1;

        RLCi[er] = j - b;
        b ^= f;
        er += f;

        ERi[j] = er;
        x1 = x0;
    }
    x0 = 0;
    f = x0 ^ x1;
    RLCi[er] = w - b;
    er += f;
    *ner = er;
}

void era_init(uint ERAi[], uf_t* EQ, uint ner) {
    for (size_t er = 1; er <= ner; er += 2)
        ERAi[er] = eq_new_label(EQ);
}

// this algorithm detects vertical neighbouring
// between RLC "segments" of current and previous row
// vertically adjacent segments are then marked as equivalent
void eq(
    uint const ERi0[],
    uint const RLCi[],
    uint const ERAi0[],
    uint ERAi[],
    uf_t* EQ,
    uint ner,
    uint W)
{
    // iterate on RLC spans
    for (size_t er = 1; er <= ner; er += 2) {
        // left end and right end of current RLC span as pixel indices
        uint j0 = RLCi[er - 1], j1 = RLCi[er];
        // uncomment for 8-connectivity
        //  j0 -= j0 > 0;
        //  j1 += j1 < W - 1;

        // leftmost and rightmost labels of the current RLC span
        // but of the previous row
        int er0 = ERi0[j0], er1 = ERi0[j1];
        // [check label parity: segments are odd]
        //  if (er0 % 2 == 0) er0++;
        //  if (er1 % 2 == 0) er1--;
            er0 += !(er0 & 1);
            er1 -= !(er1 & 1);

        // if any segments above
        if (er1 >= er0) {
            uint ea = ERAi0[er0],
                a = eq_find_root(EQ, ea);

            // missing "step 2" in the paper
            for (int erk = er0 + 2; erk <= er1; erk += 2) {
                uint eak = ERAi0[erk],
                    ak = eq_find_root(EQ, eak);

                // [min extraction and propagation]
                if (a < ak) {
                    eq_update_table(EQ, ak, a);
                }
                else if (a > ak) {
                    eq_update_table(EQ, a, ak);
                    a = ak;
                }
            }
            ERAi[er] = a; // [the global min]
        }
        else {
            ERAi[er] = eq_new_label(EQ);
        }
    }
}

void label(
    uint Ai[], uint w,
    uint const ERAi[], uint const ERi[],
    uf_t* EQ)
{
    for (size_t j = 0; j < w; ++j) {
        Ai[j] = eq_get_label(EQ, ERAi[ERi[j]]);
    }
}

size_t lsl(
    Xt const in[], uint out[],
    lsl_t const* a)
{
    size_t W = a->W, H = a->H, na = 0;//, nrlc = 0;
    uint
        ERAw    = W + 1,
        RLCw    = (W + 1) & ~1,
        *ER     = a->ER,
        *ERA    = a->ERA,
        *RLC    = a->RLC,
        *NER    = a->NER;
    //  *NRLC   = a->NRLC;
    uf_t EQ;
        EQ.EQ   = a->EQ;
        EQ.n    = 0;

    for (size_t y = 0; y < H; ++y) {
        seg(in + W * y, W,
            ER + W * y,
            RLC + RLCw * y, //RLC + nrlc,
            NER + y);
    //  NRLC[y] = nrlc;
    //  nrlc += NER[y];
    }
    memset(ERA, 0, ERAw * H * sizeof(uint));
    era_init(ERA, &EQ, NER[0]);
    for (size_t y = 1; y < H; ++y) {
        eq(ER + W * (y - 1),
            RLC + RLCw * y, //RLC + NRLC[y],
            ERA + ERAw * (y - 1),
            ERA + ERAw * y,
            &EQ, NER[y], W);
    }
    na = eq_flatten(&EQ);
    for (size_t y = 0; y < H; ++y) {
        label(out + W * y, W,
            ERA + ERAw * y,
            ER + W * y, &EQ);
    }
    return na;
}

size_t lsl_tables_size(size_t W, size_t H) {
    return
        W * H +             // ER
        (W + 1) * H +       // ERA
        ((W + 1) & ~1) * H +    // RLC
        H +                 // NER
    //  H +                 // NRLC
        W * H;              // EQ
}

lsl_t lsl_tables_from_ptr(void* memory, size_t W, size_t H) {
    lsl_t m;
        m.ER    = (uint*) memory;
        m.ERA   = m.ER      + W * H;
        m.RLC   = m.ERA     + (W + 1) * H;
        m.NER   = m.RLC     + ((W + 1) & ~1) * H;
    //  m.NRLC  = m.NER     + H;
        m.EQ    = m.NER     + H; //m.EQ     = m.NRLC    + H;
        m.W     = W;
        m.H     = H;
    return m;
}



uint eq_new_label(uf_t* uf) {
    uf->EQ[uf->n] = uf->n;
    return uf->n++;
}

uint eq_get_label(uf_t const* uf, uint i) {
    return uf->EQ[i];
}

uint eq_find_root(uf_t const* uf, uint er) {
    while (uf->EQ[er] < er) {
        er = uf->EQ[er];
    }
    return er;
}

void eq_update_table(uf_t* uf, uint e, uint r) {
    uf->EQ[e] = r;
}

uint eq_flatten(uf_t* uf) {
    uint k = 1, *EQ = uf->EQ;
    for (uint i = 1; i < uf->n; ++i) {
        EQ[i] = (EQ[i] != i)
            ?   EQ[EQ[i]]
            :   k++;
    }
    return k;
}