diff --git a/src/roi/Roi.ts b/src/roi/Roi.ts index 9eb8b3842..8ec3ca6ee 100644 --- a/src/roi/Roi.ts +++ b/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; } } diff --git a/src/roi/__tests__/ellipse.test.ts b/src/roi/__tests__/ellipse.test.ts new file mode 100644 index 000000000..0ab0a142a --- /dev/null +++ b/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, + }); +}); diff --git a/src/roi/properties/getEllipse.ts b/src/roi/properties/getEllipse.ts new file mode 100644 index 000000000..a3caaf964 --- /dev/null +++ b/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, + }; +}