feat: add cropRectangle (#439)

Closes: #436

Author: stropitek

src/Image.ts

@@ -37,12 +37,12 @@ import {
   GradientFilterOptions,
   hypotenuse,
   HypotenuseOptions,
+  increaseContrast,
+  IncreaseContrastOptions,
   invert,
   InvertOptions,
   level,
   LevelOptions,
-  increaseContrast,
-  IncreaseContrastOptions,
   medianFilter,
   MedianFilterOptions,
   pixelate,
@@ -86,14 +86,16 @@ import {
   copyTo,
   CopyToOptions,
   crop,
+  cropAlpha,
   CropAlphaOptions,
   CropOptions,
+  cropRectangle,
+  CropRectangleOptions,
   grey,
   paintMaskOnImage,
   PaintMaskOnImageOptions,
   split,
 } from './operations';
-import { cropAlpha } from './operations/cropAlpha';
 import { colorModels, ImageColorModel } from './utils/constants/colorModels';
 import { getMinMax } from './utils/getMinMax';
 import { validateChannel, validateValue } from './utils/validators/validators';
@@ -842,6 +844,17 @@ export class Image {
     return crop(this, options);
   }
 
+  /**
+   * Crop an oriented rectangle from the image.
+   * If the rectangle's length or width are not an integers, its dimension is expanded in both directions such as the length and width are integers.
+   * @param points - The points of the rectangle. Points must be circling around the rectangle (clockwise or anti-clockwise)
+   * @param options - Crop options, see {@link CropRectangleOptions}
+   * @returns The cropped image. The orientation of the image is the one closest to the rectangle passed as input.
+   */
+  public cropRectangle(points: Point[], options?: CropRectangleOptions) {
+    return cropRectangle(this, points, options);
+  }
+
   /**
    * Crops the image based on the alpha channel
    * This removes lines and columns where the alpha channel is lower than a threshold value.

src/geometry/transform.ts

@@ -2,7 +2,7 @@ import { inverse, Matrix } from 'ml-matrix';
 
 import { Image } from '../Image';
 import { getClamp } from '../utils/clamp';
-import { getBorderInterpolation, BorderType } from '../utils/interpolateBorder';
+import { BorderType, getBorderInterpolation } from '../utils/interpolateBorder';
 import {
   getInterpolationFunction,
   InterpolationType,

src/maskAnalysis/utils/getAngle.ts

@@ -2,12 +2,12 @@ import { difference, normalize, Point } from '../../utils/geometry/points';
 
 /**
  * The angle in radians of a vector relatively to the x axis.
- * The angle is positive in the cloxkwise direction.
+ * The angle is positive in the clockwise direction.
  * This is an optimized version because it assumes that one of
  * the points is on the line y = 0.
  * @param p1 - First point.
  * @param p2 - Second point.
- * @returns Rotation angle in radians to make the line horizontal.
+ * @returns Rotation angle in radians to make the line horizontal. -π <= angle <= π.
  */
 export function getAngle(p1: Point, p2: Point): number {
   const diff = difference(p2, p1);

src/operations/__tests__/cropRectangle.test.ts

@@ -0,0 +1,145 @@
+import { Image } from '../../Image';
+import { rotatePoint } from '../../point/operations';
+import { Point } from '../../utils/geometry/points';
+import { assert } from '../../utils/validators/assert';
+
+test('straight rectangle top left', () => {
+  const image = testUtils.createGreyImage([
+    [1, 2, 3],
+    [4, 5, 6],
+    [7, 8, 9],
+  ]);
+  const points = testUtils.createPoints([0, 0], [2, 0], [2, 2], [0, 2]);
+
+  const result = image.cropRectangle(points, { interpolationType: 'nearest' });
+  expect(result).toMatchImage(
+    testUtils.createGreyImage([
+      [1, 2],
+      [4, 5],
+    ]),
+  );
+});
+
+test('straight rectangle bottom right', () => {
+  const image = testUtils.createGreyImage([
+    [1, 2, 3],
+    [4, 5, 6],
+    [7, 8, 9],
+  ]);
+  const points = testUtils.createPoints([1, 3], [3, 3], [3, 1], [1, 1]);
+
+  const result = image.cropRectangle(points);
+  expect(result).toMatchImage(
+    testUtils.createGreyImage([
+      [5, 6],
+      [8, 9],
+    ]),
+  );
+});
+
+test('vertical rectangle with small angle', () => {
+  const image = testUtils.createGreyImage([
+    [1, 2, 3],
+    [4, 5, 6],
+    [7, 8, 9],
+  ]);
+  const points = testUtils
+    .createPoints([1.5, 0], [2.5, 0], [2.5, 3], [1.5, 3])
+    .map((p) => rotatePoint(p, { row: 1.5, column: 1.5 }, 0.1));
+
+  const expected = testUtils.createGreyImage([[3], [6], [8]]);
+  expectCropRectangleToMatch({
+    image,
+    expected,
+    points,
+  });
+});
+
+test('horizontal rectangle with small angle', () => {
+  const image = testUtils.createGreyImage([
+    [1, 2, 3],
+    [4, 5, 6],
+    [7, 8, 9],
+  ]);
+  const points = testUtils
+    .createPoints([0, 1.5], [3, 1.5], [3, 2.5], [0, 2.5])
+    .map((p) => rotatePoint(p, { row: 1.5, column: 1.5 }, 0.1));
+
+  const expected = testUtils.createGreyImage([[4, 5, 9]]);
+  expectCropRectangleToMatch({
+    image,
+    expected,
+    points,
+  });
+});
+
+test('diagonal rectangle oriented slightly > 45 degrees clockwise', () => {
+  const image = testUtils.createGreyImage([
+    [1, 2, 3],
+    [4, 5, 6],
+    [7, 8, 9],
+  ]);
+  const points = testUtils
+    .createPoints([1.5, 0], [2.5, 0], [2.5, 3], [1.5, 3])
+    .map((p) => rotatePoint(p, { row: 1.5, column: 1.5 }, Math.PI / 4 + 0.01));
+
+  // Resulting image is horizontal because the rectangle is closer to the horizontal axis than the vertical one
+  const expected = testUtils.createGreyImage([[0, 8, 6]]);
+  expectCropRectangleToMatch({
+    image,
+    expected,
+    points,
+  });
+});
+
+test('diagonal rectangle oriented slightly below 45 degrees', () => {
+  const image = testUtils.createGreyImage([
+    [1, 2, 3],
+    [4, 5, 6],
+    [7, 8, 9],
+  ]);
+
+  const points = testUtils
+    .createPoints([1.5, 0], [2.5, 0], [2.5, 3], [1.5, 3])
+    .map((p) => rotatePoint(p, { row: 1.5, column: 1.5 }, Math.PI / 4 - 0.01));
+
+  // Resulting image is vertical because the rectangle is closer to the vertical axis than the horizontal one
+  const expected = testUtils.createGreyImage([[0], [6], [8]]);
+  expectCropRectangleToMatch({
+    image,
+    expected,
+    points,
+  });
+});
+
+function expectCropRectangleToMatch(options: {
+  image: Image;
+  expected: Image;
+  points: Point[];
+}) {
+  const { image, expected, points } = options;
+  assert(options.points.length === 4, 'Expected to receive 4 points');
+  const variants: Point[][] = [];
+  const pointsReversed = points.slice().reverse();
+  for (let i = 0; i < 4; i++) {
+    const variant: Point[] = [];
+    variants.push(variant);
+    for (let j = 0; j < 4; j++) {
+      variant.push(points[(i + j) % 4]);
+    }
+  }
+  for (let i = 0; i < 4; i++) {
+    const variant: Point[] = [];
+    variants.push(variant);
+    for (let j = 0; j < 4; j++) {
+      variant.push(pointsReversed[(i + j) % 4]);
+    }
+  }
+
+  for (const points of variants) {
+    const result = image.cropRectangle(points, {
+      interpolationType: 'nearest',
+    });
+    expect(result).toMatchImage(expected);
+  }
+}

src/operations/cropRectangle.ts

@@ -0,0 +1,111 @@
+import { Image } from '../Image';
+import { transform, TransformOptions } from '../geometry';
+import { getAngle } from '../maskAnalysis/utils/getAngle';
+import { rotatePoint } from '../point/operations';
+import { Point } from '../utils/geometry/points';
+
+export type CropRectangleOptions = Omit<
+  TransformOptions,
+  'width' | 'height' | 'inverse' | 'fullImage'
+>;
+
+/**
+ * Crop an oriented rectangle from an image.
+ * If the rectangle's length or width are not an integers, its dimension is expanded in both directions such as the length and width are integers.
+ * @param image - The input image
+ * @param points - The points of the rectangle. Points must be circling around the rectangle (clockwise or anti-clockwise). The validity of the points passed is assumed and not checked.
+ * @param options - Crop options, see {@link CropRectangleOptions}
+ * @returns The cropped image. The orientation of the image is the one closest to the rectangle passed as input.
+ */
+export function cropRectangle(
+  image: Image,
+  points: Point[],
+  options?: CropRectangleOptions,
+): Image {
+  if (points.length !== 4) {
+    throw new Error('The points array must contain 4 points');
+  }
+
+  // get the smallest possible angle which puts the rectangle in an upright position
+  const angle = getSmallestAngle(points);
+
+  const center: Point = {
+    row: (points[0].row + points[2].row) / 2,
+    column: (points[0].column + points[2].column) / 2,
+  };
+
+  // Rotated points form an upright rectangle
+  const rotatedPoints = points.map((p) => rotatePoint(p, center, angle));
+  const [p1, p2, p3] = rotatedPoints;
+
+  const originalWidth = Math.max(
+    Math.abs(p1.column - p2.column),
+    Math.abs(p2.column - p3.column),
+  );
+  const originalHeight = Math.max(
+    Math.abs(p1.row - p2.row),
+    Math.abs(p2.row - p3.row),
+  );
+
+  // Deal with numerical imprecision when the rectangle actually had a whole number width or height
+  const width = Math.min(
+    Math.ceil(originalWidth),
+    Math.ceil(originalWidth - 1e-10),
+  );
+  const height = Math.min(
+    Math.ceil(originalHeight),
+    Math.ceil(originalHeight - 1e-10),
+  );
+
+  // Top left position of the upright rectangle after normalization of width and height
+  const expandedTopLeft = {
+    row:
+      Math.min(...rotatedPoints.map((p) => p.row)) -
+      (height - originalHeight) / 2,
+    column:
+      Math.min(...rotatedPoints.map((p) => p.column)) -
+      (width - originalWidth) / 2,
+  };
+
+  const translation = rotatePoint(expandedTopLeft, center, -angle);
+
+  const angleCos = Math.cos(-angle);
+  const angleSin = Math.sin(-angle);
+  const matrix = [
+    [angleCos, -angleSin, translation.column],
+    [angleSin, angleCos, translation.row],
+  ];
+
+  return transform(image, matrix, {
+    inverse: true,
+    width,
+    height,
+    ...options,
+  });
+}
+
+/**
+ * Get the smallest angle to put the rectangle in an upright position
+ * @param points - 2 points forming a line
+ * @returns The angle in radians
+ */
+function getSmallestAngle(points: Point[]): number {
+  // Angle respective to horizontal, -π/2 and π/2
+  let angleHorizontal = -getAngle(points[1], points[0]);
+
+  if (angleHorizontal > Math.PI / 2) {
+    angleHorizontal -= Math.PI;
+  } else if (angleHorizontal < -Math.PI / 2) {
+    angleHorizontal += Math.PI;
+  }
+
+  // angle is between -π/4 and π/4
+  let angle = angleHorizontal;
+  if (Math.abs(angleHorizontal) > Math.PI / 4) {
+    angle =
+      angleHorizontal > 0
+        ? -Math.PI / 2 + angleHorizontal
+        : Math.PI / 2 + angleHorizontal;
+  }
+  return angle;
+}

src/operations/index.ts

@@ -7,6 +7,8 @@ export * from './threshold';
 export * from './grey';
 export * from './copyTo';
 export * from './crop';
+export * from './cropAlpha';
+export * from './cropRectangle';
 export * from './operations.types';
 export * from './paintMaskOnImage';
 export * from './paintMaskOnMask';

src/point/operations.ts

@@ -0,0 +1,29 @@
+import { Point } from '../utils/geometry/points';
+
+/**
+ * Rotate a point around a center by a given angle.
+ * @param point - The point to rotate
+ * @param rotationCenter - The center of rotation
+ * @param angle - The angle of rotation in radians
+ * @returns The rotated point
+ */
+export function rotatePoint(
+  point: Point,
+  rotationCenter: Point,
+  angle: number,
+): Point {
+  const angleCos = Math.cos(angle);
+  const angleSin = Math.sin(angle);
+
+  const column =
+    point.column * angleCos -
+    point.row * angleSin +
+    (1 - angleCos) * rotationCenter.column +
+    rotationCenter.row * angleSin;
+  const row =
+    point.column * angleSin +
+    point.row * angleCos +
+    (1 - angleCos) * rotationCenter.row -
+    rotationCenter.column * angleSin;
+  return { column, row };
+}

src/stack/compute/meanImage.ts

@@ -1,6 +1,5 @@
 import { Image } from '../../Image';
 import { Stack } from '../../Stack';
-
 import { checkProcessable } from '../utils/checkProcessable';
 
 /**

src/stack/compute/medianImage.ts

@@ -3,7 +3,6 @@ import quickMedian from 'median-quickselect';
 
 import { Image } from '../../Image';
 import { Stack } from '../../Stack';
-
 import { checkProcessable } from '../utils/checkProcessable';
 
 /**

src/stack/compute/minImage.ts

@@ -1,6 +1,5 @@
 import { Image } from '../../Image';
 import { Stack } from '../../Stack';
-
 import { checkProcessable } from '../utils/checkProcessable';
 
 /**