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);
        }
    }
}