246 new roi property ellipse (#309)

* wip: add ellipse as a new feature

* chore: get ellipse to scale to match ROI's surface

* refactor: change surface name to lower case

* test: fix testing case result

* refactor: change names of the variables and return types

* chore: remove unnecessary packages

* refactor: change the way of notation of power

* test: fix the testing case to match previous commits

* chore: redo the whole implementation of an ellipse property due to a bug

* test: change test case results due to refactoring

* chore: remove commented import

* fix: add a missing bracket to border interface

* fix: remove ts-expect errors and fix object property names

* chore: remove unused dependency

* test: add testing cases for better code coverage

* chore: convert angle in ellipse from rad to degrees

* test: convert angles to degrees in testing cases

* chore: remove unused ml-array-variance

* chore: move function getEllipse to a separate folder

* refactor: refactor conditions for eigenvalues

* docs: add some comments about function and its parameters

* chore: scale the ellipse internally instead of calculating the ellipse, comparing it and scaling it

* test: fix minorAxis expecting  result

* chore: remove surface as parameter and nbSD since no longer used

* test: add testing case for 1 1 1 mask

Author: EscapedGibbon

src/roi/Roi.ts

@@ -12,6 +12,7 @@ import { Point } from '../utils/geometry/points';
 import { RoiMap } from './RoiMapManager';
 import { getBorderPoints } from './getBorderPoints';
 import { getMask, GetMaskOptions } from './getMask';
+import { Ellipse, getEllipse } from './properties/getEllipse';
 
 interface Border {
   connectedID: number; // refers to the roiID of the contiguous ROI
@@ -36,6 +37,7 @@ interface Computed {
   fillRatio: number;
   internalIDs: number[];
   feret: Feret;
+  ellipse: Ellipse;
   centroid: Point;
 }
 export class Roi {
@@ -392,6 +394,13 @@ export class Roi {
     });
   }
 
+  get ellipse(): Ellipse {
+    return this.#getComputed('ellipse', () => {
+      const ellipse = getEllipse(this);
+      return ellipse;
+    });
+  }
+
   /**
    * Number of holes in the ROI and their total surface.
    * Used to calculate fillRatio.
@@ -608,7 +617,7 @@ export class Roi {
    * @param x
    */
   computeIndex(y: number, x: number): number {
-    const roiMap = this.getMap();
+    const roiMap = this.map;
     return (y + this.origin.row) * roiMap.width + x + this.origin.column;
   }
 }

src/roi/__tests__/ellipse.test.ts

@@ -0,0 +1,160 @@
+import { fromMask } from '..';
+
+test('ellipse on a small figure 3x3', () => {
+  const mask = testUtils.createMask([
+    [0, 1, 0],
+    [0, 1, 0],
+    [0, 1, 0],
+  ]);
+  const roiMapManager = fromMask(mask);
+  const rois = roiMapManager.getRois();
+  const result = rois[0].ellipse;
+
+  expect(result).toBeDeepCloseTo({
+    center: { column: 1, row: 1 },
+    majorAxis: {
+      points: [
+        { column: NaN, row: Infinity },
+        { column: NaN, row: -Infinity },
+      ],
+      length: Infinity,
+      angle: NaN,
+    },
+    minorAxis: {
+      points: [
+        { column: NaN, row: NaN },
+        { column: NaN, row: NaN },
+      ],
+      length: NaN,
+      angle: NaN,
+    },
+    surface: NaN,
+  });
+});
+test('ellipse on a small figure 3x3', () => {
+  const mask = testUtils.createMask([
+    [1, 1, 0],
+    [0, 1, 0],
+    [0, 0, 0],
+  ]);
+  const roiMapManager = fromMask(mask);
+
+  const rois = roiMapManager.getRois();
+  const result = rois[0].ellipse;
+
+  expect(result).toBeDeepCloseTo({
+    center: { column: 0.6666666666666666, row: 0.3333333333333333 },
+    majorAxis: {
+      points: [
+        { column: 1.4978340587735344, row: 1.1645007254402011 },
+        { column: -0.1645007254402011, row: -0.4978340587735344 },
+      ],
+      length: 2.3508963970396173,
+      angle: -135,
+    },
+    minorAxis: {
+      points: [
+        { column: 1.146541384241206, row: -0.14654138424120605 },
+        { column: 0.18679194909212726, row: 0.8132080509078727 },
+      ],
+      length: 1.6247924149339357,
+      angle: 135,
+    },
+    surface: 3,
+  });
+});
+
+test('ellipse on 3x3 cross', () => {
+  const mask = testUtils.createMask([
+    [0, 1, 0],
+    [1, 1, 1],
+    [0, 1, 0],
+  ]);
+  const roiMapManager = fromMask(mask);
+
+  const rois = roiMapManager.getRois();
+  const result = rois[0].ellipse;
+  expect(result).toBeDeepCloseTo({
+    center: { column: 1, row: 1 },
+    majorAxis: {
+      points: [
+        { column: 1, row: 2.7841241161527712 },
+        { column: 1, row: -0.7841241161527714 },
+      ],
+      length: 3.5682482323055424,
+      angle: -90,
+    },
+    minorAxis: {
+      points: [
+        { column: 2.7841241161527712, row: 1 },
+        { column: -0.7841241161527714, row: 1 },
+      ],
+      length: 1.7841241161527712,
+      angle: 180,
+    },
+    surface: 5,
+  });
+});
+test('ellipse on slightly changed 3x3 cross', () => {
+  const mask = testUtils.createMask([
+    [1, 1, 0],
+    [1, 1, 1],
+    [0, 1, 0],
+  ]);
+  const roiMapManager = fromMask(mask);
+
+  const rois = roiMapManager.getRois();
+  const result = rois[0].ellipse;
+  expect(result).toBeDeepCloseTo({
+    center: { column: 0.8333333333333334, row: 0.8333333333333334 },
+    majorAxis: {
+      points: [
+        { column: 2.175183801009782, row: 2.175183801009782 },
+        { column: -0.5085171343431153, row: -0.5085171343431155 },
+      ],
+      length: 3.7953262601294284,
+      angle: -135,
+    },
+    minorAxis: {
+      points: [
+        { column: -0.15768891509232064, row: 1.8243555817589874 },
+        { column: 1.8243555817589874, row: -0.15768891509232064 },
+      ],
+      length: 2.01285390103734,
+      angle: -45,
+    },
+    surface: 6.000000000000002,
+  });
+});
+test('ellipse on 4x4 ROI', () => {
+  const mask = testUtils.createMask([
+    [0, 0, 1, 1],
+    [0, 0, 1, 0],
+    [0, 1, 1, 1],
+    [1, 1, 1, 0],
+  ]);
+  const roiMapManager = fromMask(mask);
+
+  const rois = roiMapManager.getRois();
+  const result = rois[0].ellipse;
+  expect(result).toBeDeepCloseTo({
+    center: { column: 1.7777777777777777, row: 1.7777777777777777 },
+    majorAxis: {
+      points: [
+        { column: 0.488918397751106, row: 3.6243064249674615 },
+        { column: 3.0666371578044496, row: -0.06875086941190611 },
+      ],
+      length: 4.503699166851579,
+      angle: -55.08532670592521,
+    },
+    minorAxis: {
+      points: [
+        { column: 0.8646151602330147, row: 1.1403989822180598 },
+        { column: 2.690940395322541, row: 2.4151565733374953 },
+      ],
+      length: 2.544387508597132,
+      angle: 34.91467329407479,
+    },
+    surface: 9.000000000000002,
+  });
+});

src/roi/properties/getEllipse.ts

@@ -0,0 +1,135 @@
+import { EigenvalueDecomposition } from 'ml-matrix';
+import { xVariance, xyCovariance } from 'ml-spectra-processing';
+
+import { FeretDiameter } from '../../maskAnalysis';
+import { getAngle } from '../../maskAnalysis/utils/getAngle';
+import { assert } from '../../utils/assert';
+import { toDegrees } from '../../utils/geometry/angles';
+import { Roi } from '../Roi';
+
+export interface Ellipse {
+  center: {
+    column: number;
+    row: number;
+  };
+  majorAxis: FeretDiameter;
+  minorAxis: FeretDiameter;
+  surface: number;
+}
+/**
+ *Calculates ellipse on around ROI
+ *
+ * @param roi - region of interest
+ * @param surface - the surface of ROI that ellipse should match
+ * @returns Ellipse
+ */
+export function getEllipse(roi: Roi): Ellipse {
+  let xCenter = roi.centroid.column;
+  let yCenter = roi.centroid.row;
+
+  let xCentered = roi.points.map((point: number[]) => point[0] - xCenter);
+  let yCentered = roi.points.map((point: number[]) => point[1] - yCenter);
+
+  let centeredXVariance = xVariance(xCentered, { unbiased: false });
+  let centeredYVariance = xVariance(yCentered, { unbiased: false });
+
+  let centeredCovariance = xyCovariance(
+    {
+      x: xCentered,
+      y: yCentered,
+    },
+    { unbiased: false },
+  );
+
+  //spectral decomposition of the sample covariance matrix
+  let sampleCovarianceMatrix = [
+    [centeredXVariance, centeredCovariance],
+    [centeredCovariance, centeredYVariance],
+  ];
+  let e = new EigenvalueDecomposition(sampleCovarianceMatrix);
+  let eigenvalues = e.realEigenvalues;
+  let vectors = e.eigenvectorMatrix;
+
+  let radiusMajor: number;
+  let radiusMinor: number;
+  let vectorMajor: number[];
+  let vectorMinor: number[];
+
+  assert(eigenvalues[0] <= eigenvalues[1]);
+  radiusMajor = Math.sqrt(eigenvalues[1]);
+  radiusMinor = Math.sqrt(eigenvalues[0]);
+  vectorMajor = vectors.getColumn(1);
+  vectorMinor = vectors.getColumn(0);
+
+  let majorAxisPoint1 = {
+    column: xCenter + radiusMajor * vectorMajor[0],
+    row: yCenter + radiusMajor * vectorMajor[1],
+  };
+  let majorAxisPoint2 = {
+    column: xCenter - radiusMajor * vectorMajor[0],
+    row: yCenter - radiusMajor * vectorMajor[1],
+  };
+  let minorAxisPoint1 = {
+    column: xCenter + radiusMinor * vectorMinor[0],
+    row: yCenter + radiusMinor * vectorMinor[1],
+  };
+  let minorAxisPoint2 = {
+    column: xCenter - radiusMinor * vectorMinor[0],
+    row: yCenter - radiusMinor * vectorMinor[1],
+  };
+
+  let majorLength = Math.sqrt(
+    (majorAxisPoint1.column - majorAxisPoint2.column) ** 2 +
+      (majorAxisPoint1.row - majorAxisPoint2.row) ** 2,
+  );
+  let minorLength = Math.sqrt(
+    (minorAxisPoint1.column - majorAxisPoint2.column) ** 2 +
+      (minorAxisPoint1.row - minorAxisPoint2.row) ** 2,
+  );
+
+  let ellipseSurface = (((minorLength / 2) * majorLength) / 2) * Math.PI;
+  if (ellipseSurface !== roi.surface) {
+    const scaleFactor = Math.sqrt(roi.surface / ellipseSurface);
+    radiusMajor *= scaleFactor;
+    radiusMinor *= scaleFactor;
+    majorAxisPoint1 = {
+      column: xCenter + radiusMajor * vectorMajor[0],
+      row: yCenter + radiusMajor * vectorMajor[1],
+    };
+    majorAxisPoint2 = {
+      column: xCenter - radiusMajor * vectorMajor[0],
+      row: yCenter - radiusMajor * vectorMajor[1],
+    };
+    minorAxisPoint1 = {
+      column: xCenter + radiusMinor * vectorMinor[0],
+      row: yCenter + radiusMinor * vectorMinor[1],
+    };
+    minorAxisPoint2 = {
+      column: xCenter - radiusMinor * vectorMinor[0],
+      row: yCenter - radiusMinor * vectorMinor[1],
+    };
+
+    majorLength *= scaleFactor;
+
+    minorLength *= scaleFactor;
+    ellipseSurface *= scaleFactor ** 2;
+  }
+
+  return {
+    center: {
+      column: xCenter,
+      row: yCenter,
+    },
+    majorAxis: {
+      points: [majorAxisPoint1, majorAxisPoint2],
+      length: majorLength,
+      angle: toDegrees(getAngle(majorAxisPoint1, majorAxisPoint2)),
+    },
+    minorAxis: {
+      points: [minorAxisPoint1, minorAxisPoint2],
+      length: minorLength,
+      angle: toDegrees(getAngle(minorAxisPoint1, minorAxisPoint2)),
+    },
+    surface: ellipseSurface,
+  };
+}