Connected-component out of area

[KnewHow]

I view your algorithms in playne_equivalence_direct_3D_clamp.cu and I implement with dx11 computer shader with same code

cbuffer Dim : register(b0)
{
    uint X;
    uint Y;
    uint Z;
    uint XYZ;
    uint PX;
    uint PY;
    uint useless1;
    uint useless2;
};

StructuredBuffer<uint> Volume: register(t0); // A volume with uint to indicate Connected-Component
RWStructuredBuffer<uint> Labels: register(u0); // A lable has same size with Volume to tag which field is connected

// ---------- Find the root of a chain ----------
inline uint find_root(uint label) {
    // Resolve Label
    uint next = Labels[label];

    // Follow chain
    while (label != next) {
        // Move to next
        label = next;
        next = Labels[label];
    }

    // Return label
    return label;
}

// ---------- Label Reduction ----------
inline uint reduction(uint label1, uint label2) {
    // Get next labels
    uint next1 = (label1 != label2) ? Labels[label1] : 0;
    uint next2 = (label1 != label2) ? Labels[label2] : 0;

    // Find label1
    while ((label1 != label2) && (label1 != next1)) {
        // Adopt label
        label1 = next1;

        // Fetch next label
        next1 = Labels[label1];
    }

    // Find label2
    while ((label1 != label2) && (label2 != next2)) {
        // Adopt label
        label2 = next2;

        // Fetch next label
        next2 = Labels[label2];
    }

    uint label3;
    // While Labels are different
    while (label1 != label2) {
        // Label 2 should be smallest
        if (label1 < label2) {
            // Swap Labels
            label1 = label1 ^ label2;
            label2 = label1 ^ label2;
            label1 = label1 ^ label2;
        }

        // AtomicMin label1 to label2
        InterlockedMin(Labels[label1], label2, label3);
        label1 = (label1 == label3) ? label2 : label3;
    }

    // Return label1
    return label1;
}

[numthreads(4, 4, 4)]
void initLabel(
    uint3 DTid : SV_DispatchThreadID) 
{
    const uint ix = DTid.x;
    const uint iy = DTid.y;
    const uint iz = DTid.z;
    if ((ix < X) && (iy < Y) && (iz < Z)) {
        const uint pzyx = Volume[iz * PY + iy * PX + ix];
        // Neighbour Connections
        const bool nzm1yx = (iz > 0) ? (pzyx == Volume[(iz - 1) * PY + iy * PX + ix]) : false;
        const bool nzym1x = (iy > 0) ? (pzyx == Volume[iz * PY + (iy - 1) * PX + ix]) : false;
        const bool nzyxm1 = (ix > 0) ? (pzyx == Volume[iz * PY + iy * PX + ix - 1]) : false;

        // Label
        uint label;

        // Initialise Label
        label = (nzyxm1) ? (iz * PY + iy * PX + ix - 1) : (iz * PY + iy * PX + ix);
        label = (nzym1x) ? (iz * PY + (iy - 1) * PX + ix) : label;
        label = (nzm1yx) ? ((iz - 1) * PY + iy * PX + ix) : label;

        // Write to Global Memory
        Labels[iz * PY + iy * PX + ix] = label;
    }
}

[numthreads(4, 4, 4)]
void resolveLabel(
    uint3 DTid : SV_DispatchThreadID) 
{
    const uint ix = DTid.x;
    const uint iy = DTid.y;
    const uint iz = DTid.z;
    // Calculate index
    const uint idx = iz * PY + iy * PX + ix;

    // Check Range
    if (idx < XYZ) {
        // Resolve Label
        Labels[idx] = find_root(Labels[idx]);
    }
}

[numthreads(4, 4, 4)]
void labelReduction(
    uint3 DTid : SV_DispatchThreadID) 
{
    // Calculate index
    const uint ix = DTid.x;
    const uint iy = DTid.y;
    const uint iz = DTid.z;

    if (ix < X && iy < Y && iz < Z) {
        const uint pzyx = Volume[iz * PY + iy * PX + ix];

        // Compare Image Values
        const bool nzm1yx = (iz > 0) ? (pzyx == Labels[(iz - 1) * PY + iy * PX + ix]) : false;
        const bool nzym1x = (iy > 0) ? (pzyx == Labels[iz * PY + (iy - 1) * PX + ix]) : false;
        const bool nzyxm1 = (ix > 0) ? (pzyx == Labels[iz * PY + iy * PX + ix - 1]) : false;

        const bool nzym1xm1 = ((iy > 0) && (ix > 0)) ? (pzyx == Labels[iz * PY + (iy - 1) * PX + ix - 1]) : false;
        const bool nzm1yxm1 = ((iz > 0) && (ix > 0)) ? (pzyx == Labels[(iz - 1) * PY + iy * PX + ix - 1]) : false;
        const bool nzm1ym1x = ((iz > 0) && (iy > 0)) ? (pzyx == Labels[(iz - 1) * PY + (iy - 1) * PX + ix]) : false;

        /*unsigned int label1 = (nzm1yx || nzym1x || nzyxm1) ? Labels[iz * PY + iy * PX + ix] : 0;
        unsigned int label2 = 0;
        if (nzm1yx) {
            unsigned int label2 = Labels[(iz - 1) * PY + iy * PX + ix];
            label1 = reduction(label1, label2);
        }
        if (nzym1x) {
            unsigned int label2 = Labels[iz * PY + (iy - 1) * PX + ix];
            label1 = reduction(label1, label2);
        }
        if (nzyxm1) {
            unsigned int label2 = Labels[iz * PY + iy * PX + ix - 1];
            label1 = reduction(label1, label2);
        }
        return;*/

        // Critical point conditions
        const bool cond_x = nzyxm1 && ((nzm1yx && !nzm1yxm1) || (nzym1x && !nzym1xm1));
        const bool cond_y = nzym1x && nzm1yx && !nzm1ym1x;

        // Get label
        uint label1 = (cond_x || cond_y) ? Labels[iz * PY + iy * PX + ix] : 0;

        // Y - Neighbour
        if (cond_y) {
            // Get neighbouring label
            uint label2 = Labels[iz * PY + (iy - 1) * PX + ix];

            // Reduction
            label1 = reduction(label1, label2);
        }

        // X - Neighbour
        if (cond_x) {
            // Get neighbouring label
            uint label2 = Labels[iz * PY + iy * PX + ix - 1];

            // Reduction
            label1 = reduction(label1, label2);
        }
    }

}

[numthreads(4, 4, 4)]
void clean(
    uint3 DTid : SV_DispatchThreadID) 
{
    const uint ix = DTid.x;
    const uint iy = DTid.y;
    const uint iz = DTid.z;
    Labels[iz * PY + iy * PX + ix] = 0;
}

Then I test it, I have a 128x128x128 cube volume, and I init it with a sphere in center, I think the connected-component is sphere volume, So I init the cube volume with 1 in shpere area, otherwise is zero.

int dimX = 128;
            int dimY = 128;
            int dimZ = 128;
            float vxSizeX = 1.0f / dimX;
            float vxSizeY = 1.0f / dimY;
            float vxSizeZ = 1.0f / dimZ;

            float originX = 0;
            float originY = 0;
            float originZ = 0;

            float centerX = dimX / 2;
            float centerY = dimY / 2;
            float centerZ = dimZ / 2;
            float radius = (dimX / 2.7f) * (dimX / 2.7f);
            byte[] data = new byte[dimX * dimY * dimZ];
            for (int k = 0; k < dimZ; k++)
            {
                for (int j = 0; j < dimY; j++)
                {
                    for (int i = 0; i < dimX; i++)
                    {
                        long index = k * dimX * dimY + j * dimX + i;
                        if ((i - centerX) * (i - centerX) + (j - centerY) * (j - centerY) + (k - centerZ) * (k - centerZ) < radius)
                        {
                            data[index] = 128;
                        }
                        else
                        {
                            data[index] = 0;
                        }

                    }
                }
            }
            float[] grid = new float[9] { originX, originY, originZ, vxSizeX, vxSizeY, vxSizeZ, dimX, dimY, dimZ };
            return new Volume(grid, data);

Then I run the algorithms, and get the Labels, then I print the lable with following code:

int counter = 0;
            StreamWriter writer = new StreamWriter("C:\\temp\\a.txt");
            for (int k = 0; k < dimZ; k++)
            {
                for (int j = 0; j < dimY; j++)
                {
                    for (int i = 0; i < dimX; i++)
                    {
                        long index = k * dimX * dimY + j * dimX + i;
                        long idx = index * 2;
                        writer.Write(labels[index] + " ");
                        ++counter;
                        if(counter % 32 == 0)
                        {
                            writer.Write("\n");
                        }
                        if (labels[index] > 0)
                        {
                            data[idx] = 1;
                            data[idx + 1] = 1;
                        } else
                        {
                            data[idx] = 0;
                            data[idx + 1] = 0;
                        }

                    }
                }
            }
            writer.Flush();
            writer.Close();

The I print it each 32 with a newline, you can download the file from this link: Lables data Then I use volume render to render it with above init: if the Labes value in current index is greater than zero, I give the voxel with two bytes in [1,1], otherwise I give the voxel with two bytes is[0,0], Then I use accumulated volume render to render it, I find some question in it. You can get the render image from this: https://drive.google.com/file/d/1KY5mQ5JcsRzA4bPxCRNh8tSPeuYPa2rG/view?usp=sharing Why some connected-component outer of the sphere, I think the connected-component should in the inner of the shpere. Is there some bug in my code? Would you help me to solve this problem?

[DanielPlayne]

Hi @KnewHow,

Thanks for you enquiry. I don't have an environment set up to test it but I suspect the problem occurs on lines 10-16 of your labelReduction function where you compute the neighbouring connections - nzm1yx, nzym1x, nzyxm1, nzym1xm1, nzm1yxm1, nzm1ym1x.

These should represent the connections between neighbouring voxels and in the original implementation it tests whether the image value pzyx is equal to the image value of various neighbours. Your code seems to compare the value of the image pzyx with the value of the neighbouring labels instead.

What this means is that you'll only be getting labels from the init_label step where the kernel looks at the neighbours and initialises the label with the best guess from the neighbours.

As another note, this implementation of component labelling will compute a label for every voxel in your dataset, it doesn't discriminate between background and foreground image data as I originally wrote it for some statistical mechanics. If you want it to ignore the background voxels (data points with a value of 0) then you would need to add a condition to the init_label and label_reduction to ignore voxels with a value of 0.

Regards

• Daniel

[KnewHow]

Thank you very much, you solve my problem. Best wish! My code has some problem in

 const bool nzm1yx = (iz > 0) ? (pzyx == Labels[(iz - 1) * PY + iy * PX + ix]) : false;
        const bool nzym1x = (iy > 0) ? (pzyx == Labels[iz * PY + (iy - 1) * PX + ix]) : false;
        const bool nzyxm1 = (ix > 0) ? (pzyx == Labels[iz * PY + iy * PX + ix - 1]) : false;

        const bool nzym1xm1 = ((iy > 0) && (ix > 0)) ? (pzyx == Labels[iz * PY + (iy - 1) * PX + ix - 1]) : false;
        const bool nzm1yxm1 = ((iz > 0) && (ix > 0)) ? (pzyx == Labels[(iz - 1) * PY + iy * PX + ix - 1]) : false;
        const bool nzm1ym1x = ((iz > 0) && (iy > 0)) ? (pzyx == Labels[(iz - 1) * PY + (iy - 1) * PX + ix]) : false;

It should be

 const bool nzm1yx = (iz > 0) ? (pzyx == Volume[(iz - 1) * PY + iy * PX + ix]) : false;
        const bool nzym1x = (iy > 0) ? (pzyx == Volume[iz * PY + (iy - 1) * PX + ix]) : false;
        const bool nzyxm1 = (ix > 0) ? (pzyx == Volume[iz * PY + iy * PX + ix - 1]) : false;

        const bool nzym1xm1 = ((iy > 0) && (ix > 0)) ? (pzyx == Volume[iz * PY + (iy - 1) * PX + ix - 1]) : false;
        const bool nzm1yxm1 = ((iz > 0) && (ix > 0)) ? (pzyx == Volume[(iz - 1) * PY + iy * PX + ix - 1]) : false;
        const bool nzm1ym1x = ((iz > 0) && (iy > 0)) ? (pzyx == Volume[(iz - 1) * PY + (iy - 1) * PX + ix]) : false;