export class BWLabeler {
// port of https://github.com/rordenlab/niimath/blob/master/src/bwlabel.c
// return voxel address given row A, column B, and slice C
idx(A, B, C, DIM) {
return C * DIM[0] * DIM[1] + B * DIM[0] + A
} // idx()
// determine if voxels below candidate voxel have already been assigned a label
check_previous_slice(bw, il, r, c, sl, dim, conn, tt, nabo, tn) {
let nr_set = 0
if (!sl) {
return 0
}
const val = bw[this.idx(r, c, sl, dim)]
if (conn >= 6) {
const idx = this.idx(r, c, sl - 1, dim)
if (val === bw[idx]) {
nabo[nr_set++] = il[idx]
}
}
if (conn >= 18) {
if (r) {
const idx = this.idx(r - 1, c, sl - 1, dim)
if (val === bw[idx]) {
nabo[nr_set++] = il[idx]
}
}
if (c) {
const idx = this.idx(r, c - 1, sl - 1, dim)
if (val === bw[idx]) {
nabo[nr_set++] = il[idx]
}
}
if (r < dim[0] - 1) {
const idx = this.idx(r + 1, c, sl - 1, dim)
if (val === bw[idx]) {
nabo[nr_set++] = il[idx]
}
}
if (c < dim[1] - 1) {
const idx = this.idx(r, c + 1, sl - 1, dim)
if (val === bw[idx]) {
nabo[nr_set++] = il[idx]
}
}
}
if (conn === 26) {
if (r && c) {
const idx = this.idx(r - 1, c - 1, sl - 1, dim)
if (val === bw[idx]) {
nabo[nr_set++] = il[idx]
}
}
if (r < dim[0] - 1 && c) {
const idx = this.idx(r + 1, c - 1, sl - 1, dim)
if (val === bw[idx]) {
nabo[nr_set++] = il[idx]
}
}
if (r && c < dim[1] - 1) {
const idx = this.idx(r - 1, c + 1, sl - 1, dim)
if (val === bw[idx]) {
nabo[nr_set++] = il[idx]
}
}
if (r < dim[0] - 1 && c < dim[1] - 1) {
const idx = this.idx(r + 1, c + 1, sl - 1, dim)
if (val === bw[idx]) {
nabo[nr_set++] = il[idx]
}
}
}
if (nr_set) {
this.fill_tratab(tt, nabo, nr_set, tn)
return nabo[0]
} else {
return 0
}
} // check_previous_slice()
// provisionally label all voxels in volume
do_initial_labelling(bw, dim, conn) {
const naboPS = new Uint32Array(32)
const tn = new Uint32Array(32)
let label = 1
const kGrowArrayBy = 8192
let ttn = kGrowArrayBy
let tt = new Uint32Array(ttn).fill(0)
const il = new Uint32Array(dim[0] * dim[1] * dim[2]).fill(0)
const nabo = new Uint32Array(27)
for (let sl = 0; sl < dim[2]; sl++) {
for (let c = 0; c < dim[1]; c++) {
for (let r = 0; r < dim[0]; r++) {
let nr_set = 0
const val = bw[this.idx(r, c, sl, dim)]
if (val === 0) {
continue
}
nabo[0] = this.check_previous_slice(bw, il, r, c, sl, dim, conn, tt, naboPS, tn)
if (nabo[0]) {
nr_set += 1
}
if (conn >= 6) {
if (r) {
const idx = this.idx(r - 1, c, sl, dim)
if (val === bw[idx]) {
nabo[nr_set++] = il[idx]
}
}
if (c) {
const idx = this.idx(r, c - 1, sl, dim)
if (val === bw[idx]) {
nabo[nr_set++] = il[idx]
}
}
}
if (conn >= 18) {
if (c && r) {
const idx = this.idx(r - 1, c - 1, sl, dim)
if (val === bw[idx]) {
nabo[nr_set++] = il[idx]
}
}
if (c && r < dim[0] - 1) {
const idx = this.idx(r + 1, c - 1, sl, dim)
if (val === bw[idx]) {
nabo[nr_set++] = il[idx]
}
}
}
if (nr_set) {
il[this.idx(r, c, sl, dim)] = nabo[0]
this.fill_tratab(tt, nabo, nr_set, tn)
} else {
il[this.idx(r, c, sl, dim)] = label
if (label >= ttn) {
ttn += kGrowArrayBy
const ext = new Uint32Array(ttn)
ext.set(tt)
tt = ext
}
tt[label - 1] = label
label++
}
}
}
}
for (let i = 0; i < label - 1; i++) {
let j = i
while (tt[j] !== j + 1) {
j = tt[j] - 1
}
tt[i] = j + 1
}
return [label - 1, tt, il]
} // do_initial_labelling()
// translation table unifies a region that has been assigned multiple classes
fill_tratab(tt, nabo, nr_set, tn) {
// let cntr = 0
//tn.fill(0)
const INT_MAX = 2147483647
let ltn = INT_MAX
for (let i = 0; i < nr_set; i++) {
let j = nabo[i]
// cntr = 0
while (tt[j - 1] !== j) {
j = tt[j - 1]
/* cntr++
if (cntr > 100) {
console.log('\nOoh no!!')
break
} */
}
tn[i] = j
ltn = Math.min(ltn, j)
}
for (let i = 0; i < nr_set; i++) {
tt[tn[i] - 1] = ltn
}
} // fill_tratab()
// remove any residual gaps so label numbers are dense rather than sparse
translate_labels(il, dim, tt, ttn) {
const nvox = dim[0] * dim[1] * dim[2]
let ml = 0
const l = new Uint32Array(nvox).fill(0)
for (let i = 0; i < ttn; i++) {
ml = Math.max(ml, tt[i])
}
const fl = new Uint32Array(ml).fill(0)
let cl = 0
for (let i = 0; i < nvox; i++) {
if (il[i]) {
if (!fl[tt[il[i] - 1] - 1]) {
cl += 1
fl[tt[il[i] - 1] - 1] = cl
}
l[i] = fl[tt[il[i] - 1] - 1]
}
}
return [cl, l]
} // translate_labels()
// retain only the largest cluster for each region
largest_original_cluster_labels(bw, cl, ls) {
const nvox = bw.length
const ls2bw = new Uint32Array(cl + 1).fill(0)
const sumls = new Uint32Array(cl + 1).fill(0)
for (let i = 0; i < nvox; i++) {
const bwVal = bw[i]
const lsVal = ls[i]
ls2bw[lsVal] = bwVal
sumls[lsVal]++
}
let mxbw = 0
for (let i = 0; i < cl + 1; i++) {
const bwVal = ls2bw[i]
mxbw = Math.max(mxbw, bwVal)
// see if this is largest cluster of this bw-value
for (let j = 0; j < cl + 1; j++) {
if (j === i) {
continue
}
if (bwVal !== ls2bw[j]) {
continue
}
if (sumls[i] < sumls[j]) {
ls2bw[i] = 0
} else if (sumls[i] === sumls[j] && i < j) {
ls2bw[i] = 0
} // ties: arbitrary winner
}
}
const vxs = new Uint32Array(nvox).fill(0)
for (let i = 0; i < nvox; i++) {
vxs[i] = ls2bw[ls[i]]
}
return [mxbw, vxs]
}
// given a 3D image, return a clustered label map
// for an explanation and optimized C code see
// https://github.com/seung-lab/connected-components-3d
bwlabel(img, dim, conn = 26, binarize = false, onlyLargestClusterPerClass = false) {
const start = Date.now()
const nvox = dim[0] * dim[1] * dim[2]
const bw = new Uint32Array(nvox).fill(0)
if (![6, 18, 26].includes(conn)) {
console.log('bwlabel: conn must be 6, 18 or 26.')
return [0, bw]
}
if (dim[0] < 2 || dim[1] < 2 || dim[2] < 1) {
console.log('bwlabel: img must be 2 or 3-dimensional')
return [0, bw]
}
if (binarize) {
for (let i = 0; i < nvox; i++) {
if (img[i] !== 0.0) {
bw[i] = 1
}
}
} else {
bw.set(img)
}
let [ttn, tt, il] = this.do_initial_labelling(bw, dim, conn)
if (tt === undefined) {
tt = new Uint32Array(0)
}
const [cl, ls] = this.translate_labels(il, dim, tt, ttn)
console.log(conn + ' neighbor clustering into ' + cl + ' regions in ' + (Date.now() - start) + 'ms')
if (onlyLargestClusterPerClass) {
const [nbw, bwMx] = this.largest_original_cluster_labels(bw, cl, ls)
return [nbw, bwMx]
}
return [cl, ls]
} // bwlabel()
}
