# Lindenberg's watershed-based blob detection for MATLAB

Just something I needed recently and haven’t found any existing standalone implementation. Algorithm description can be found on Wikipedia. Usage is quite simple:

% Where I is 1-, 2- or 3-dimensional image
[~, index] = sort(reshape(I, numel(I), 1), 1, ‘descend’);
blobs = ls_blob_detect(I, index);

```/* Lindenberg's watershed-based blob detection algorithm implementation for MATLAB
Author: Stanislaw Adaszewski, 2012 */

#include <mex.h>

static void ind2sub(int nx, int ny, int nz, int idx, int *x, int *y, int *z) {
int r1 = idx % (nx * ny);
*z = idx / (nx * ny);
*y = r1 / nx;
*x = r1 % nx;
}

static int sub2ind(int nx, int ny, int nz, int x, int y, int z) {
if (x < 0 || x >= nx || y < 0 || y >= ny || z < 0 || z >= nz)
return -1;
else
return x + nx * (y + ny * z);
}

static const int ofs[][3] = {
{-1, -1, -1}, {-1,  0, -1}, {-1, 1, -1}, {0, -1, -1}, {0, 0, -1}, {0, 1, -1}, {1, -1, -1}, {1, 0, -1}, {1, 1, -1},
{-1, -1,  0}, {-1,  0,  0}, {-1, 1,  0}, {0, -1,  0},
// {0, 0, 0},
{ 0,  1,  0}, { 1, -1,  0}, { 1, 0,  0}, {1,  1,  0},
{-1, -1,  1}, {-1,  0,  1}, {-1, 1,  1}, {0, -1,  1}, {0, 0,  1}, {0, 1,  1}, {1, -1,  1}, {1, 0,  1}, {1, 1,  1}
};

void mexFunction (int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
if (nlhs != 1 || nrhs != 2 || mxGetNumberOfDimensions(prhs[0]) > 3 || !mxIsDouble(prhs[0]) || mxGetClassID(prhs[1]) != mxINT32_CLASS || mxGetNumberOfElements(prhs[1]) != mxGetNumberOfElements(prhs[0])) {
mexErrMsgTxt("Usage: R=ls_blob_detect(I, index)nI must be 1-, 2- or 3-dimensional imagenindex must be list of voxel offsets (0-based) in decreasing brightness order");
return;
} else {
int N = mxGetNumberOfElements(prhs[0]);
double *data1 = mxGetPr(prhs[0]);
int n_blobs = 0;
const mwSize *siz = mxGetDimensions(prhs[0]);
int ndim = mxGetNumberOfDimensions(prhs[0]);
mxArray *R = mxCreateNumericArray(ndim, siz, mxINT32_CLASS, mxREAL);
int nx, ny, nz;
int *map = (int*) mxGetData(R);
int *index = (int*) mxGetData(prhs[1]); // malloc(sizeof(int) * N);
int i;

if (!R) {
mexErrMsgTxt("Out of memory.");
return;
}

plhs[0] = R;

for (i = 0; i < N; i++)
map[i] = 0;

nx = siz[0];
ny = siz[1];
if (ndim > 2)
nz = siz[2];
else
nz = 1;

for (i = 0; i < N; i++) {
int x, y, z;
int *neighbors = (int*) calloc(sizeof(int), n_blobs+1);
int j;
int blob, blob_cnt = 0;
if (!neighbors) {
mexErrMsgTxt("Out of memory.");
return;
}

if (data1[index[i]] == 0)
continue; // Ignore 0 voxels

ind2sub(nx, ny, nz, index[i], &x, &y, &z);

for (j = 0; j < 26; j++) {
int idx_neighbor = sub2ind(nx, ny, nz, x + ofs[j][0], y + ofs[j][1], z + ofs[j][2]);

if (idx_neighbor == -1)
continue; // Ignore voxels outside of volume boundaries

if (data1[idx_neighbor] > data1[index[i]])
neighbors[map[idx_neighbor]] += 1;
}

for (j = 0; j < n_blobs; j++)
if (neighbors[1+j] > 0) {
blob_cnt += 1;
blob = j;
}

if (neighbors[0]>0 || blob_cnt>1) {
map[index[i]] = 0; // background
} else if (blob_cnt > 0) {
map[index[i]] = 1 + blob; // existing blob
} else {
n_blobs += 1;
map[index[i]] = n_blobs; // new blob
}

free(neighbors);
}
}
}```

• Cless says:

My input is a 2-D grayscalce image. I follow the usage precisely and keep getting this error below.

So what’s wrong?

Thank you for sharing !

“Error using ls_blob_detect
Usage: R=ls_blob_detect(I, index)
I must be 1-, 2- or 3-dimensional image
index must be list of voxel offsets (0-based) in decreasing brightness order”