feat: add watershed filter (#379)

* feat: add watershed filter

* feat: add threshold computation to watershed filter

* test: add a testing case for threshold

* fix: update package.json

* fix: fix eslint errors

* fix: fix eslint errors in testing file

* test: add test for coverage

* fix: fix the problem with return values and improve extrema check

* test: add testing cases

* fix: fix eslint error

* fix: remove unnecessary points check

* fix: fix RoiMapManager computation

* fix: fix eslint error

* refactor: add kind:'bw' to get all the existing rois

* feat: add a helper function to get Int16Arrays for easier presentation and  debugging

* docs: add documentation to getInt16Array function

* test: refactor output arrays for better presentation and simpler debugging and change test names

* refactor: remove 'kind' option and adjust the code

* chore: move watershed to roi folder

* chore: add watershed export to index.ts

* chore: fix eslint errors

Author: EscapedGibbon

package.json

@@ -51,6 +51,7 @@
     "fast-png": "^6.2.0",
     "image-type": "^4.1.0",
     "jpeg-js": "^0.4.4",
+    "js-priority-queue": "^0.1.5",
     "median-quickselect": "^1.0.1",
     "ml-affine-transform": "^1.0.3",
     "ml-convolution": "^2.0.0",
@@ -66,6 +67,7 @@
     "@microsoft/api-extractor": "^7.36.3",
     "@tailwindcss/forms": "^0.5.4",
     "@types/jest": "^29.5.3",
+    "@types/js-priority-queue": "^0.0.6",
     "@types/jest-image-snapshot": "^6.2.0",
     "@types/node": "^20.4.8",
     "@types/picomatch": "^2.3.0",

src/roi/RoiMapManager.ts

@@ -81,7 +81,7 @@ export class RoiMapManager implements RoiManager {
   }
 
   public getRoiById(roiID: number): Roi {
-    const rois = this.getRois();
+    const rois = this.getRois({ kind: 'bw' });
     const foundRoi = rois.find((roi) => roi.id === roiID);
     if (!foundRoi) {
       throw new Error(`invalid ID: ${roiID}`);

src/roi/__tests__/waterShed.test.ts

@@ -0,0 +1,197 @@
+import { computeThreshold } from '../..';
+import { createGreyImage, getInt16Array } from '../../../test/testUtils';
+import { waterShed } from '../waterShed';
+
+describe('Test WaterShed Roi generation', () => {
+  it('test 1,basic test without parameters/options', () => {
+    const image = createGreyImage([
+      [3, 3, 3, 3, 3],
+      [3, 2, 2, 2, 3],
+      [3, 2, 1, 2, 3],
+      [3, 2, 2, 2, 3],
+      [3, 3, 3, 3, 3],
+    ]);
+    const roiMapManager = waterShed(image, {});
+    const resultArray = testUtils.getInt16Array(`
+      -1, -1, -1, -1, -1,
+      -1, -1, -1, -1, -1,
+      -1, -1, -1, -1, -1,
+      -1, -1, -1, -1, -1,
+      -1, -1, -1, -1, -1,`);
+    expect(roiMapManager).toEqual({
+      map: {
+        data: resultArray,
+        nbPositive: 0,
+        nbNegative: 1,
+        width: 5,
+        height: 5,
+      },
+      whiteRois: [],
+      blackRois: [],
+    });
+  });
+
+  it('test 2, waterShed for a grey image', () => {
+    const image = createGreyImage([
+      [3, 3, 3, 3, 3, 3, 3, 3, 4, 4],
+      [3, 3, 2, 2, 2, 3, 3, 3, 4, 4],
+      [4, 3, 2, 1, 2, 2, 3, 3, 4, 4],
+      [4, 3, 2, 2, 2, 2, 3, 3, 3, 4],
+      [4, 4, 4, 3, 2, 3, 2, 3, 3, 4],
+      [4, 4, 4, 3, 3, 3, 3, 1, 3, 3],
+      [4, 3, 3, 3, 3, 3, 2, 2, 2, 3],
+      [4, 4, 3, 3, 3, 3, 2, 2, 2, 2],
+      [4, 4, 4, 4, 3, 2, 2, 2, 2, 3],
+      [4, 4, 4, 4, 3, 3, 3, 3, 2, 3],
+    ]);
+
+    const roiMapManager = waterShed(image, { threshold: 2 / 255 });
+
+    const resultArray = getInt16Array(`
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+      0, 0,-2,-2,-2, 0, 0, 0, 0, 0,
+      0, 0,-2,-2,-2,-2, 0, 0, 0, 0,
+      0, 0,-2,-2,-2,-2, 0, 0, 0, 0,
+      0, 0, 0, 0,-2, 0, 0, 0, 0, 0,
+      0, 0, 0, 0, 0, 0, 0,-1, 0, 0,
+      0, 0, 0, 0, 0, 0,-1,-1,-1, 0,
+      0, 0, 0, 0, 0, 0,-1,-1,-1,-1,
+      0, 0, 0, 0, 0,-1,-1,-1,-1, 0,
+      0, 0, 0, 0, 0, 0, 0, 0,-1, 0,`);
+    expect(roiMapManager).toEqual({
+      map: {
+        data: resultArray,
+        nbPositive: 0,
+        nbNegative: 2,
+        width: 10,
+        height: 10,
+      },
+      whiteRois: [],
+      blackRois: [],
+    });
+    const rois = roiMapManager.getRois({ kind: 'bw' });
+    expect(rois[0].origin).toEqual({ column: 5, row: 5 });
+    expect(rois[1].origin).toEqual({ column: 2, row: 1 });
+  });
+
+  it('test 3, with threshold option', () => {
+    const image = createGreyImage([
+      [1, 1, 1, 1, 1],
+      [1, 2, 2, 2, 1],
+      [1, 2, 3, 2, 1],
+      [1, 2, 2, 2, 1],
+      [1, 1, 1, 1, 1],
+    ]);
+    const invertedImage = image.invert();
+    const roiMapManager = waterShed(invertedImage, {
+      threshold: 253 / image.maxValue,
+    });
+    const resultArray = getInt16Array(`
+      0, 0, 0, 0, 0,
+      0,-1,-1,-1, 0, 
+      0,-1,-1,-1, 0, 
+      0,-1,-1,-1, 0,
+      0, 0, 0, 0, 0,
+    `);
+    expect(roiMapManager).toEqual({
+      map: {
+        data: resultArray,
+        nbPositive: 0,
+        nbNegative: 1,
+        width: 5,
+        height: 5,
+      },
+      whiteRois: [],
+      blackRois: [],
+    });
+    const roi1 = roiMapManager.getRoiById(-1);
+
+    expect(roi1.surface).toEqual(9);
+    expect(roi1.width).toEqual(3);
+  });
+  it('test 4, waterShed through threshold value', () => {
+    const image = createGreyImage([
+      [3, 3, 3, 3, 3, 3, 3, 3, 4, 4],
+      [3, 3, 2, 2, 2, 3, 3, 3, 4, 4],
+      [4, 3, 2, 1, 2, 2, 3, 3, 4, 4],
+      [4, 3, 2, 2, 2, 2, 3, 3, 3, 4],
+      [4, 4, 4, 3, 2, 3, 3, 3, 3, 4],
+      [4, 4, 4, 3, 3, 3, 3, 3, 3, 3],
+      [4, 3, 3, 3, 3, 3, 2, 2, 2, 3],
+      [4, 4, 3, 3, 3, 3, 2, 1, 2, 2],
+      [4, 4, 4, 4, 3, 2, 2, 2, 2, 3],
+      [4, 4, 4, 4, 3, 3, 3, 3, 2, 3],
+    ]);
+    const threshold = computeThreshold(image, 'otsu');
+
+    const roiMapManager = waterShed(image, {
+      threshold: threshold / image.maxValue,
+    });
+    const resultArray = getInt16Array(`
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+      0, 0,-2,-2,-2, 0, 0, 0, 0, 0, 
+      0, 0,-2,-2,-2,-2, 0, 0, 0, 0,
+      0, 0,-2,-2,-2,-2, 0, 0, 0, 0,
+      0, 0, 0, 0,-2, 0, 0, 0, 0, 0,
+      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
+      0, 0, 0, 0, 0, 0,-1,-1,-1, 0,
+      0, 0, 0, 0, 0, 0,-1,-1,-1,-1,
+      0, 0, 0, 0, 0,-1,-1,-1,-1, 0,
+      0, 0, 0, 0, 0, 0, 0, 0,-1, 0,
+    `);
+    expect(roiMapManager).toEqual({
+      map: {
+        data: resultArray,
+        nbPositive: 0,
+        nbNegative: 2,
+        width: 10,
+        height: 10,
+      },
+      whiteRois: [],
+      blackRois: [],
+    });
+  });
+  it('test 5, waterShed through threshold mask and with inverted image', () => {
+    const image = createGreyImage([
+      [3, 3, 3, 3, 3, 3, 3, 3, 4, 4],
+      [3, 3, 2, 2, 2, 3, 3, 8, 4, 4],
+      [4, 3, 2, 1, 2, 2, 3, 3, 4, 4],
+      [4, 3, 2, 2, 6, 2, 3, 3, 3, 4],
+      [4, 4, 4, 3, 2, 3, 3, 3, 3, 4],
+      [4, 4, 4, 3, 3, 3, 3, 3, 3, 3],
+      [4, 3, 3, 3, 3, 6, 2, 2, 2, 3],
+      [4, 4, 3, 3, 3, 3, 2, 1, 2, 2],
+      [4, 4, 4, 4, 3, 2, 2, 2, 2, 3],
+      [4, 4, 4, 4, 3, 3, 3, 3, 9, 3],
+    ]);
+
+    const mask = image.threshold({ algorithm: 'otsu' });
+    const roiMapManager = waterShed(image, { mask });
+    const resultArray = getInt16Array(`
+      -2, -2, -2, -2, -2, -2, -2, -2, -2, -1,
+      -2, -2, -2, -2, -2, -2, -2,  0, -1, -1,
+      -2, -2, -2, -2, -2, -2, -2, -2, -1, -1, 
+      -2, -2, -2, -2,  0, -2, -2, -1, -1, -1,
+      -2, -2, -2, -2, -1, -2, -1, -1, -1, -1,
+      -2, -1, -1, -1, -1, -1, -1, -1, -1, -1,
+      -1, -1, -1, -1, -1,  0, -1, -1, -1, -1,
+      -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
+      -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
+      -1, -1, -1, -1, -1, -1, -1, -1,  0, -1,
+    `);
+    expect(roiMapManager).toEqual({
+      map: {
+        data: resultArray,
+        nbPositive: 0,
+        nbNegative: 2,
+        width: 10,
+        height: 10,
+      },
+      whiteRois: [],
+      blackRois: [],
+    });
+    const rois = roiMapManager.getRois({ kind: 'bw' });
+    expect(rois[0].origin).toEqual({ column: 0, row: 0 });
+    expect(rois[0].surface).toEqual(60);
+  });
+});

src/roi/computeRois.ts

@@ -37,7 +37,9 @@ export function computeRois(roiMapManager: RoiMapManager): void {
   for (let row = 0; row < map.height; row++) {
     for (let column = 0; column < map.width; column++) {
       const currentIndex = roiMapManager.getMapValue(column, row);
-
+      if (currentIndex === 0) {
+        continue;
+      }
       let currentRoi;
       if (currentIndex < 0) {
         currentRoi = blacks[-currentIndex - 1];

src/roi/index.ts

@@ -6,3 +6,4 @@ export * from './getMask';
 export * from './getRois';
 export * from './Roi';
 export * from './RoiMapManager';
+export * from './waterShed';

src/roi/waterShed.ts

@@ -0,0 +1,146 @@
+import PriorityQueue from 'js-priority-queue';
+
+import { RoiMapManager } from '..';
+import { Image } from '../Image';
+import { Mask } from '../Mask';
+import getExtrema from '../compute/getExtrema';
+import { Point } from '../geometry';
+import { filterPoints } from '../utils/geometry/filterPoints';
+import checkProcessable from '../utils/validators/checkProcessable';
+/**
+ * Point interface that is used in the queue data structure.
+ */
+interface PointWithIntensity {
+  /**
+   * @param row - Row of a point.
+   */
+  row: number;
+  /**
+   * @param column - Column of a point.
+   */
+  column: number;
+  /**
+   * @param intensity - Value of a point.
+   */
+  intensity: number;
+}
+
+interface WaterShedOptions {
+  /**
+   * @param points - Points which should be filled by watershed filter.
+   * @default - minimum points from getExtrema() function.
+   */
+  points?: Point[];
+  /**
+   * @param mask - A binary image, the same size as the image. The algorithm will fill only if the current pixel in the binary mask is not null.
+   * @default undefined
+   */
+  mask?: Mask;
+  /**
+   * @param threshold - Limit of filling. Maximum value that pixel can have.
+   * @default 1
+   */
+  threshold?: number;
+}
+/**
+ * This method allows to create a ROIMap using the water shed algorithm. By default this algorithm
+ * will fill the holes and therefore the lowest value of the image (black zones).
+ * If no points are given, the function will look for all the minimal points.
+ * If no mask is given the algorithm will completely fill the image.
+ * Please take care about the value that has be in the mask ! In order to be coherent with the expected mask,
+ * meaning that if it is a dark zone, the mask will be dark the normal behavior to fill a zone
+ * is that the mask pixel is clear (value of 0) !
+ * If you are looking for 'maxima' the image must be inverted before applying the algorithm
+ * @param image - Image that the filter will be applied to.
+ * @param options - WaterShedOptions
+ * @returns RoiMapManager
+ */
+export function waterShed(
+  image: Image,
+  options: WaterShedOptions,
+): RoiMapManager {
+  let { points } = options;
+  const { mask, threshold = 1 } = options;
+  const currentImage = image;
+  checkProcessable(image, {
+    bitDepth: [8, 16],
+    components: 1,
+  });
+
+  const fillMaxValue = threshold * image.maxValue;
+
+  // WaterShed is done from points in the image. We can either specify those points in options,
+  // or it is gonna take the minimum locals of the image by default.
+  if (!points) {
+    points = getExtrema(image, {
+      kind: 'minimum',
+      mask,
+    });
+    points = filterPoints(points, image, { kind: 'minimum' });
+  }
+
+  const maskExpectedValue = 0;
+
+  const data = new Int16Array(currentImage.size);
+  const width = currentImage.width;
+  const height = currentImage.height;
+  const toProcess = new PriorityQueue({
+    comparator: (a: PointWithIntensity, b: PointWithIntensity) =>
+      a.intensity - b.intensity,
+    strategy: PriorityQueue.BinaryHeapStrategy,
+  });
+  for (let i = 0; i < points.length; i++) {
+    const index = points[i].column + points[i].row * width;
+    data[index] = -i - 1;
+    const intensity = currentImage.getValueByIndex(index, 0);
+    if (intensity <= fillMaxValue) {
+      toProcess.queue({
+        column: points[i].column,
+        row: points[i].row,
+        intensity,
+      });
+    }
+  }
+  const dxs = [+1, 0, -1, 0, +1, +1, -1, -1];
+  const dys = [0, +1, 0, -1, +1, -1, +1, -1];
+  // Then we iterate through each points
+
+  while (toProcess.length > 0) {
+    const currentPoint = toProcess.dequeue();
+    const currentValueIndex = currentPoint.column + currentPoint.row * width;
+    for (let dir = 0; dir < 4; dir++) {
+      const newX = currentPoint.column + dxs[dir];
+      const newY = currentPoint.row + dys[dir];
+      if (newX >= 0 && newY >= 0 && newX < width && newY < height) {
+        const currentNeighbourIndex = newX + newY * width;
+        if (
+          !mask ||
+          mask.getBitByIndex(currentNeighbourIndex) === maskExpectedValue
+        ) {
+          const intensity = currentImage.getValueByIndex(
+            currentNeighbourIndex,
+            0,
+          );
+          if (intensity <= fillMaxValue && data[currentNeighbourIndex] === 0) {
+            data[currentNeighbourIndex] = data[currentValueIndex];
+            toProcess.queue({
+              column: currentPoint.column + dxs[dir],
+              row: currentPoint.row + dys[dir],
+              intensity,
+            });
+          }
+        }
+      }
+    }
+  }
+  const nbNegative = points.length;
+  const nbPositive = 0;
+
+  return new RoiMapManager({
+    data,
+    nbPositive,
+    nbNegative,
+    width: image.width,
+    height: image.height,
+  });
+}