Downloads: 1
Diff of
/bwlabels.js
[000000]
..
[b86468]
Switch to unified view
| a | b/bwlabels.js | ||
|---|---|---|---|
| 1 | export class BWLabeler { |
||
| 2 | // port of https://github.com/rordenlab/niimath/blob/master/src/bwlabel.c |
||
| 3 | // return voxel address given row A, column B, and slice C |
||
| 4 | idx(A, B, C, DIM) { |
||
| 5 | return C * DIM[0] * DIM[1] + B * DIM[0] + A |
||
| 6 | } // idx() |
||
| 7 | |||
| 8 | // determine if voxels below candidate voxel have already been assigned a label |
||
| 9 | check_previous_slice(bw, il, r, c, sl, dim, conn, tt, nabo, tn) { |
||
| 10 | let nr_set = 0 |
||
| 11 | if (!sl) { |
||
| 12 | return 0 |
||
| 13 | } |
||
| 14 | const val = bw[this.idx(r, c, sl, dim)] |
||
| 15 | if (conn >= 6) { |
||
| 16 | const idx = this.idx(r, c, sl - 1, dim) |
||
| 17 | if (val === bw[idx]) { |
||
| 18 | nabo[nr_set++] = il[idx] |
||
| 19 | } |
||
| 20 | } |
||
| 21 | if (conn >= 18) { |
||
| 22 | if (r) { |
||
| 23 | const idx = this.idx(r - 1, c, sl - 1, dim) |
||
| 24 | if (val === bw[idx]) { |
||
| 25 | nabo[nr_set++] = il[idx] |
||
| 26 | } |
||
| 27 | } |
||
| 28 | if (c) { |
||
| 29 | const idx = this.idx(r, c - 1, sl - 1, dim) |
||
| 30 | if (val === bw[idx]) { |
||
| 31 | nabo[nr_set++] = il[idx] |
||
| 32 | } |
||
| 33 | } |
||
| 34 | if (r < dim[0] - 1) { |
||
| 35 | const idx = this.idx(r + 1, c, sl - 1, dim) |
||
| 36 | if (val === bw[idx]) { |
||
| 37 | nabo[nr_set++] = il[idx] |
||
| 38 | } |
||
| 39 | } |
||
| 40 | if (c < dim[1] - 1) { |
||
| 41 | const idx = this.idx(r, c + 1, sl - 1, dim) |
||
| 42 | if (val === bw[idx]) { |
||
| 43 | nabo[nr_set++] = il[idx] |
||
| 44 | } |
||
| 45 | } |
||
| 46 | } |
||
| 47 | if (conn === 26) { |
||
| 48 | if (r && c) { |
||
| 49 | const idx = this.idx(r - 1, c - 1, sl - 1, dim) |
||
| 50 | if (val === bw[idx]) { |
||
| 51 | nabo[nr_set++] = il[idx] |
||
| 52 | } |
||
| 53 | } |
||
| 54 | if (r < dim[0] - 1 && c) { |
||
| 55 | const idx = this.idx(r + 1, c - 1, sl - 1, dim) |
||
| 56 | if (val === bw[idx]) { |
||
| 57 | nabo[nr_set++] = il[idx] |
||
| 58 | } |
||
| 59 | } |
||
| 60 | if (r && c < dim[1] - 1) { |
||
| 61 | const idx = this.idx(r - 1, c + 1, sl - 1, dim) |
||
| 62 | if (val === bw[idx]) { |
||
| 63 | nabo[nr_set++] = il[idx] |
||
| 64 | } |
||
| 65 | } |
||
| 66 | if (r < dim[0] - 1 && c < dim[1] - 1) { |
||
| 67 | const idx = this.idx(r + 1, c + 1, sl - 1, dim) |
||
| 68 | if (val === bw[idx]) { |
||
| 69 | nabo[nr_set++] = il[idx] |
||
| 70 | } |
||
| 71 | } |
||
| 72 | } |
||
| 73 | if (nr_set) { |
||
| 74 | this.fill_tratab(tt, nabo, nr_set, tn) |
||
| 75 | return nabo[0] |
||
| 76 | } else { |
||
| 77 | return 0 |
||
| 78 | } |
||
| 79 | } // check_previous_slice() |
||
| 80 | |||
| 81 | // provisionally label all voxels in volume |
||
| 82 | do_initial_labelling(bw, dim, conn) { |
||
| 83 | const naboPS = new Uint32Array(32) |
||
| 84 | const tn = new Uint32Array(32) |
||
| 85 | let label = 1 |
||
| 86 | const kGrowArrayBy = 8192 |
||
| 87 | let ttn = kGrowArrayBy |
||
| 88 | let tt = new Uint32Array(ttn).fill(0) |
||
| 89 | const il = new Uint32Array(dim[0] * dim[1] * dim[2]).fill(0) |
||
| 90 | const nabo = new Uint32Array(27) |
||
| 91 | for (let sl = 0; sl < dim[2]; sl++) { |
||
| 92 | for (let c = 0; c < dim[1]; c++) { |
||
| 93 | for (let r = 0; r < dim[0]; r++) { |
||
| 94 | let nr_set = 0 |
||
| 95 | const val = bw[this.idx(r, c, sl, dim)] |
||
| 96 | if (val === 0) { |
||
| 97 | continue |
||
| 98 | } |
||
| 99 | nabo[0] = this.check_previous_slice(bw, il, r, c, sl, dim, conn, tt, naboPS, tn) |
||
| 100 | if (nabo[0]) { |
||
| 101 | nr_set += 1 |
||
| 102 | } |
||
| 103 | if (conn >= 6) { |
||
| 104 | if (r) { |
||
| 105 | const idx = this.idx(r - 1, c, sl, dim) |
||
| 106 | if (val === bw[idx]) { |
||
| 107 | nabo[nr_set++] = il[idx] |
||
| 108 | } |
||
| 109 | } |
||
| 110 | if (c) { |
||
| 111 | const idx = this.idx(r, c - 1, sl, dim) |
||
| 112 | if (val === bw[idx]) { |
||
| 113 | nabo[nr_set++] = il[idx] |
||
| 114 | } |
||
| 115 | } |
||
| 116 | } |
||
| 117 | if (conn >= 18) { |
||
| 118 | if (c && r) { |
||
| 119 | const idx = this.idx(r - 1, c - 1, sl, dim) |
||
| 120 | if (val === bw[idx]) { |
||
| 121 | nabo[nr_set++] = il[idx] |
||
| 122 | } |
||
| 123 | } |
||
| 124 | if (c && r < dim[0] - 1) { |
||
| 125 | const idx = this.idx(r + 1, c - 1, sl, dim) |
||
| 126 | if (val === bw[idx]) { |
||
| 127 | nabo[nr_set++] = il[idx] |
||
| 128 | } |
||
| 129 | } |
||
| 130 | } |
||
| 131 | if (nr_set) { |
||
| 132 | il[this.idx(r, c, sl, dim)] = nabo[0] |
||
| 133 | this.fill_tratab(tt, nabo, nr_set, tn) |
||
| 134 | } else { |
||
| 135 | il[this.idx(r, c, sl, dim)] = label |
||
| 136 | if (label >= ttn) { |
||
| 137 | ttn += kGrowArrayBy |
||
| 138 | const ext = new Uint32Array(ttn) |
||
| 139 | ext.set(tt) |
||
| 140 | tt = ext |
||
| 141 | } |
||
| 142 | tt[label - 1] = label |
||
| 143 | label++ |
||
| 144 | } |
||
| 145 | } |
||
| 146 | } |
||
| 147 | } |
||
| 148 | for (let i = 0; i < label - 1; i++) { |
||
| 149 | let j = i |
||
| 150 | while (tt[j] !== j + 1) { |
||
| 151 | j = tt[j] - 1 |
||
| 152 | } |
||
| 153 | tt[i] = j + 1 |
||
| 154 | } |
||
| 155 | return [label - 1, tt, il] |
||
| 156 | } // do_initial_labelling() |
||
| 157 | |||
| 158 | // translation table unifies a region that has been assigned multiple classes |
||
| 159 | fill_tratab(tt, nabo, nr_set, tn) { |
||
| 160 | // let cntr = 0 |
||
| 161 | //tn.fill(0) |
||
| 162 | const INT_MAX = 2147483647 |
||
| 163 | let ltn = INT_MAX |
||
| 164 | for (let i = 0; i < nr_set; i++) { |
||
| 165 | let j = nabo[i] |
||
| 166 | // cntr = 0 |
||
| 167 | while (tt[j - 1] !== j) { |
||
| 168 | j = tt[j - 1] |
||
| 169 | /* cntr++ |
||
| 170 | if (cntr > 100) { |
||
| 171 | console.log('\nOoh no!!') |
||
| 172 | break |
||
| 173 | } */ |
||
| 174 | } |
||
| 175 | tn[i] = j |
||
| 176 | ltn = Math.min(ltn, j) |
||
| 177 | } |
||
| 178 | for (let i = 0; i < nr_set; i++) { |
||
| 179 | tt[tn[i] - 1] = ltn |
||
| 180 | } |
||
| 181 | } // fill_tratab() |
||
| 182 | |||
| 183 | // remove any residual gaps so label numbers are dense rather than sparse |
||
| 184 | translate_labels(il, dim, tt, ttn) { |
||
| 185 | const nvox = dim[0] * dim[1] * dim[2] |
||
| 186 | let ml = 0 |
||
| 187 | const l = new Uint32Array(nvox).fill(0) |
||
| 188 | for (let i = 0; i < ttn; i++) { |
||
| 189 | ml = Math.max(ml, tt[i]) |
||
| 190 | } |
||
| 191 | const fl = new Uint32Array(ml).fill(0) |
||
| 192 | let cl = 0 |
||
| 193 | for (let i = 0; i < nvox; i++) { |
||
| 194 | if (il[i]) { |
||
| 195 | if (!fl[tt[il[i] - 1] - 1]) { |
||
| 196 | cl += 1 |
||
| 197 | fl[tt[il[i] - 1] - 1] = cl |
||
| 198 | } |
||
| 199 | l[i] = fl[tt[il[i] - 1] - 1] |
||
| 200 | } |
||
| 201 | } |
||
| 202 | return [cl, l] |
||
| 203 | } // translate_labels() |
||
| 204 | |||
| 205 | // retain only the largest cluster for each region |
||
| 206 | largest_original_cluster_labels(bw, cl, ls) { |
||
| 207 | const nvox = bw.length |
||
| 208 | const ls2bw = new Uint32Array(cl + 1).fill(0) |
||
| 209 | const sumls = new Uint32Array(cl + 1).fill(0) |
||
| 210 | for (let i = 0; i < nvox; i++) { |
||
| 211 | const bwVal = bw[i] |
||
| 212 | const lsVal = ls[i] |
||
| 213 | ls2bw[lsVal] = bwVal |
||
| 214 | sumls[lsVal]++ |
||
| 215 | } |
||
| 216 | let mxbw = 0 |
||
| 217 | for (let i = 0; i < cl + 1; i++) { |
||
| 218 | const bwVal = ls2bw[i] |
||
| 219 | mxbw = Math.max(mxbw, bwVal) |
||
| 220 | // see if this is largest cluster of this bw-value |
||
| 221 | for (let j = 0; j < cl + 1; j++) { |
||
| 222 | if (j === i) { |
||
| 223 | continue |
||
| 224 | } |
||
| 225 | if (bwVal !== ls2bw[j]) { |
||
| 226 | continue |
||
| 227 | } |
||
| 228 | if (sumls[i] < sumls[j]) { |
||
| 229 | ls2bw[i] = 0 |
||
| 230 | } else if (sumls[i] === sumls[j] && i < j) { |
||
| 231 | ls2bw[i] = 0 |
||
| 232 | } // ties: arbitrary winner |
||
| 233 | } |
||
| 234 | } |
||
| 235 | const vxs = new Uint32Array(nvox).fill(0) |
||
| 236 | for (let i = 0; i < nvox; i++) { |
||
| 237 | vxs[i] = ls2bw[ls[i]] |
||
| 238 | } |
||
| 239 | return [mxbw, vxs] |
||
| 240 | } |
||
| 241 | |||
| 242 | // given a 3D image, return a clustered label map |
||
| 243 | // for an explanation and optimized C code see |
||
| 244 | // https://github.com/seung-lab/connected-components-3d |
||
| 245 | bwlabel(img, dim, conn = 26, binarize = false, onlyLargestClusterPerClass = false) { |
||
| 246 | const start = Date.now() |
||
| 247 | const nvox = dim[0] * dim[1] * dim[2] |
||
| 248 | const bw = new Uint32Array(nvox).fill(0) |
||
| 249 | if (![6, 18, 26].includes(conn)) { |
||
| 250 | console.log('bwlabel: conn must be 6, 18 or 26.') |
||
| 251 | return [0, bw] |
||
| 252 | } |
||
| 253 | if (dim[0] < 2 || dim[1] < 2 || dim[2] < 1) { |
||
| 254 | console.log('bwlabel: img must be 2 or 3-dimensional') |
||
| 255 | return [0, bw] |
||
| 256 | } |
||
| 257 | if (binarize) { |
||
| 258 | for (let i = 0; i < nvox; i++) { |
||
| 259 | if (img[i] !== 0.0) { |
||
| 260 | bw[i] = 1 |
||
| 261 | } |
||
| 262 | } |
||
| 263 | } else { |
||
| 264 | bw.set(img) |
||
| 265 | } |
||
| 266 | let [ttn, tt, il] = this.do_initial_labelling(bw, dim, conn) |
||
| 267 | if (tt === undefined) { |
||
| 268 | tt = new Uint32Array(0) |
||
| 269 | } |
||
| 270 | const [cl, ls] = this.translate_labels(il, dim, tt, ttn) |
||
| 271 | console.log(conn + ' neighbor clustering into ' + cl + ' regions in ' + (Date.now() - start) + 'ms') |
||
| 272 | if (onlyLargestClusterPerClass) { |
||
| 273 | const [nbw, bwMx] = this.largest_original_cluster_labels(bw, cl, ls) |
||
| 274 | return [nbw, bwMx] |
||
| 275 | } |
||
| 276 | return [cl, ls] |
||
| 277 | } // bwlabel() |
||
| 278 | } |