/*
    * Licensed to the Apache Software Foundation (ASF) under one or more
    * contributor license agreements.  See the NOTICE file distributed with
    * this work for additional information regarding copyright ownership.
    * The ASF licenses this file to You under the Apache License, Version 2.0
    * (the "License"); you may not use this file except in compliance with
    * the License.  You may obtain a copy of the License at
    *
    *      http://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  package org.apache.commons.rng.sampling.shape;
  
  import java.util.Arrays;
  import org.apache.commons.math3.stat.inference.ChiSquareTest;
  import org.apache.commons.rng.UniformRandomProvider;
  import org.apache.commons.rng.sampling.RandomAssert;
  import org.junit.jupiter.api.Assertions;
  import org.junit.jupiter.api.Test;
  
  /**
   * Test for {@link TriangleSampler}.
   */
  class TriangleSamplerTest {
      // Precomputed 3D and 4D rotation matrices for forward and backward transforms
      // created using the matlab function by jodag:
      // https://stackoverflow.com/questions/50337642/how-to-calculate-a-rotation-matrix-in-n-dimensions-given-the-point-to-rotate-an
  
      /** 3D rotation around the vector [1; 2; 3] by pi/4. */
      private static final double[][] F3 = {
              {0.728027725387508, -0.525104821111919, 0.440727305612110},
              {0.608788597915763, 0.790790557990391, -0.063456571298848},
              {-0.315201640406345, 0.314507901710379, 0.895395278995195}};
      /** 3D rotation around the vector [1; 2; 3] by -pi/4. */
      private static final double[][] R3 = {
              {0.728027725387508, 0.608788597915762, -0.315201640406344},
              {-0.525104821111919, 0.790790557990391, 0.314507901710379},
              {0.440727305612110, -0.063456571298848, 0.895395278995195}};
      /** 4D rotation around the orthogonal vectors [1 0; 0 1; 1 0; 0 1] by pi/4. */
      private static final double[][] F4 = {
              {0.853553390593274, -0.353553390593274, 0.146446609406726, 0.353553390593274},
              {0.353553390593274, 0.853553390593274, -0.353553390593274, 0.146446609406726},
              {0.146446609406726, 0.353553390593274, 0.853553390593274, -0.353553390593274},
              {-0.353553390593274, 0.146446609406726, 0.353553390593274, 0.853553390593274}};
      /** 4D rotation around the orthogonal vectors [1 0; 0 1; 1 0; 0 1] by -pi/4. */
      private static final double[][] R4 = {
              {0.853553390593274, 0.353553390593274, 0.146446609406726, -0.353553390593274},
              {-0.353553390593274, 0.853553390593274, 0.353553390593274, 0.146446609406726},
              {0.146446609406726, -0.353553390593274, 0.853553390593274, 0.353553390593274},
              {0.353553390593274, 0.146446609406726, -0.353553390593274, 0.853553390593274}};
  
      /**
       * Test the sampling assumptions used to transform coordinates outside the triangle
       * back inside the triangle.
       */
      @Test
      void testSamplingAssumptions() {
          // The separation between the 2^53 dyadic rationals in the interval [0, 1)
          final double delta = 0x1.0p-53;
          double s = 0.5;
          double t = 0.5 + delta;
          // This value cannot be exactly represented and is rounded
          final double spt = s + t;
          // Test that (1 - (1-s) - (1-t)) is not equal to (s + t - 1).
          // This is due to the rounding to store s + t as a double.
          final double expected = 1 - (1 - s) - (1 - t);
          Assertions.assertNotEquals(expected, spt - 1);
          Assertions.assertNotEquals(expected, s + t - 1);
          // For any uniform deviate u in [0, 1], u - 1 is exact, thus s - 1 is exact
          // and s - 1 + t is exact.
          Assertions.assertEquals(expected, s - 1 + t);
  
          // Test that a(1 - s - t) + sb + tc does not overflow is s+t = 1
          final double max = Double.MAX_VALUE;
          s -= delta;
          final UniformRandomProvider rng = RandomAssert.createRNG();
          for (int n = 0; n < 100; n++) {
              Assertions.assertNotEquals(Double.POSITIVE_INFINITY, (1 - s - t) * max + s * max + t * max);
              s = rng.nextDouble();
              t = 1.0 - s;
          }
      }
  
      /**
       * Test an unsupported dimension.
       */
      @Test
      void testInvalidDimensionThrows() {
          final UniformRandomProvider rng = RandomAssert.seededRNG();
          Assertions.assertThrows(IllegalArgumentException.class,
              () -> TriangleSampler.of(rng, new double[1], new double[1], new double[1]));
      }
  
      /**
      * Test a dimension mismatch between vertices.
      */
     @Test
     void testDimensionMismatchThrows() {
         final UniformRandomProvider rng = RandomAssert.seededRNG();
         final double[] c2 = new double[2];
         final double[] c3 = new double[3];
         for (double[][] c : new double[][][] {
             {c2, c2, c3},
             {c2, c3, c2},
             {c3, c2, c2},
             {c2, c3, c3},
             {c3, c3, c2},
             {c3, c2, c3},
         }) {
             Assertions.assertThrows(IllegalArgumentException.class,
                 () -> TriangleSampler.of(rng, c[0], c[1], c[2]),
                 () -> String.format("Did not detect dimension mismatch: %d,%d,%d",
                     c[0].length, c[1].length, c[2].length));
         }
     }
 
     /**
      * Test non-finite vertices.
      */
     @Test
     void testNonFiniteVertexCoordinates() {
         final UniformRandomProvider rng = RandomAssert.seededRNG();
         // A valid triangle
         final double[][] c = new double[][] {
             {0, 0, 1}, {2, 1, 0}, {-1, 2, 3}
         };
         Assertions.assertNotNull(TriangleSampler.of(rng, c[0],  c[1],  c[2]));
         final double[] bad = {Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY, Double.NaN};
         for (int i = 0; i < c.length; i++) {
             final int ii = i;
             for (int j = 0; j < c[0].length; j++) {
                 final int jj = j;
                 for (final double d : bad) {
                     final double value = c[i][j];
                     c[i][j] = d;
                     Assertions.assertThrows(IllegalArgumentException.class,
                         () -> TriangleSampler.of(rng, c[0], c[1], c[2]),
                         () -> String.format("Did not detect non-finite coordinate: %d,%d = %s", ii, jj, d));
                     c[i][j] = value;
                 }
             }
         }
     }
 
     /**
      * Test a triangle with coordinates that are separated by more than
      * {@link Double#MAX_VALUE} in two dimensions.
      */
     @Test
     void testExtremeValueCoordinates2D() {
         testExtremeValueCoordinates(2);
     }
 
     /**
      * Test a triangle with coordinates that are separated by more than
      * {@link Double#MAX_VALUE} in three dimensions.
      */
     @Test
     void testExtremeValueCoordinates3D() {
         testExtremeValueCoordinates(3);
     }
 
     /**
      * Test a triangle with coordinates that are separated by more than
      * {@link Double#MAX_VALUE} in four dimensions.
      */
     @Test
     void testExtremeValueCoordinates4D() {
         testExtremeValueCoordinates(4);
     }
 
     /**
      * Test a triangle with coordinates that are separated by more than
      * {@link Double#MAX_VALUE}.
      *
      * @param dimension the dimension
      */
     private static void testExtremeValueCoordinates(int dimension) {
         final double[][] c1 = new double[3][dimension];
         final double[][] c2 = new double[3][dimension];
         // Create a valid triangle that can be scaled
         c1[0][0] = -1;
         c1[0][1] = -1;
         c1[1][0] = 1;
         c1[1][1] = 1;
         c1[2][0] = 1;
         c1[2][1] = -1;
         // Extremely large value for scaling. Use a power of 2 for exact scaling.
         final double scale = 0x1.0p1023;
         for (int i = 0; i < 3; i++) {
             // Fill the remaining dimensions with 1 or -1
             for (int j = 2; j < dimension; j++) {
                 c1[i][j] = 1 - (j & 0x2);
             }
             // Scale the second triangle
             for (int j = 0; j < dimension; j++) {
                 c2[i][j] = c1[i][j] * scale;
             }
         }
         // Show the triangle is too big to compute vectors between points.
         Assertions.assertEquals(Double.POSITIVE_INFINITY, c2[2][0] - c2[0][0],
             "Expect vector c - a to be infinite in the x dimension");
         Assertions.assertEquals(Double.NEGATIVE_INFINITY, c2[2][1] - c2[1][1],
             "Expect vector c - b to be infinite in the y dimension");
 
         final TriangleSampler sampler1 = TriangleSampler.of(
             RandomAssert.seededRNG(), c1[0], c1[1], c1[2]);
         final TriangleSampler sampler2 = TriangleSampler.of(
             RandomAssert.seededRNG(), c2[0], c2[1], c2[2]);
 
         for (int n = 0; n < 10; n++) {
             final double[] a = sampler1.sample();
             final double[] b = sampler2.sample();
             for (int i = 0; i < a.length; i++) {
                 a[i] *= scale;
             }
             Assertions.assertArrayEquals(a, b);
         }
     }
 
     /**
      * Test the distribution of points in two dimensions.
      */
     @Test
     void testDistribution2D() {
         testDistributionND(2);
     }
 
     /**
      * Test the distribution of points in three dimensions.
      */
     @Test
     void testDistribution3D() {
         testDistributionND(3);
     }
 
     /**
      * Test the distribution of points in four dimensions.
      */
     @Test
     void testDistribution4D() {
         testDistributionND(4);
     }
 
     /**
      * Test the distribution of points in N dimensions. The output coordinates
      * should be uniform in the triangle.
      *
      * <p>ND is supported by transforming 2D triangles to ND, sampling in ND then
      * transforming back to 2D. The transformed results should have a zero
      * coordinate for indices above 1. Supports 2 to 4 dimensions.
      *
      * @param dimension the dimension in the range [2, 4]
      */
     private static void testDistributionND(int dimension) {
         // Create 3 triangles that form a 3x4 rectangle:
         // b-----c
         // |\   /|
         // | \ / |
         // |  e  |
         // |   \ |
         // |    \|
         // a-----d
         //
         // Sample from them proportional to the area:
         // adb = 4 * 3 / 2 = 6
         // bce = 3 * 2 / 2 = 3
         // cde = 4 * 1.5 / 2 = 3
         // The samples in the ratio 2:1:1 should be uniform within the rectangle.
         final double[] a = {0, 0};
         final double[] b = {0, 4};
         final double[] c = {3, 4};
         final double[] d = {3, 0};
         final double[] e = {1.5, 2};
 
         // Assign bins
         final int bins = 20;
         final int samplesPerBin = 10;
         // Scale factors to assign x,y to a bin
         final double sx = bins / 3.0;
         final double sy = bins / 4.0;
 
         // Expect a uniform distribution (this is rescaled by the ChiSquareTest)
         final double[] expected = new double[bins * bins];
         Arrays.fill(expected, 1);
 
         // Support ND by applying a rotation transform to the 2D triangles to generate
         // n-coordinates. The samples are transformed back and should be uniform in 2D.
         Transform forward;
         Transform reverse;
         if (dimension == 4) {
             forward = new ForwardTransform(F4);
             reverse = new ReverseTransform(R4);
         } else if (dimension == 3) {
             forward = new ForwardTransform(F3);
             reverse = new ReverseTransform(R3);
         } else if (dimension == 2) {
             forward = reverse = new Transform2Dto2D();
         } else {
             throw new AssertionError("Unsupported dimension: " + dimension);
         }
 
         // Increase the loops to verify robustness
         final UniformRandomProvider rng = RandomAssert.createRNG();
         final TriangleSampler sampler1 = TriangleSampler.of(rng, forward.apply(a), forward.apply(d), forward.apply(b));
         final TriangleSampler sampler2 = TriangleSampler.of(rng, forward.apply(b), forward.apply(c), forward.apply(e));
         final TriangleSampler sampler3 = TriangleSampler.of(rng, forward.apply(c), forward.apply(d), forward.apply(e));
         final Triangle triangle1 = new Triangle(a, d, b);
         final Triangle triangle2 = new Triangle(b, c, e);
         final Triangle triangle3 = new Triangle(c, d, e);
         final int samples = expected.length * samplesPerBin;
         for (int n = 0; n < 1; n++) {
             // Assign each coordinate to a region inside the combined rectangle
             final long[] observed = new long[expected.length];
             // Sample according to the area of each triangle (ratio 2:1:1)
             for (int i = 0; i < samples; i += 4) {
                 addObservation(reverse.apply(sampler1.sample()), observed, bins, sx, sy, triangle1);
                 addObservation(reverse.apply(sampler1.sample()), observed, bins, sx, sy, triangle1);
                 addObservation(reverse.apply(sampler2.sample()), observed, bins, sx, sy, triangle2);
                 addObservation(reverse.apply(sampler3.sample()), observed, bins, sx, sy, triangle3);
             }
             final double p = new ChiSquareTest().chiSquareTest(expected, observed);
             Assertions.assertFalse(p < 0.001, () -> "p-value too small: " + p);
         }
     }
 
     /**
      * Adds the observation. Coordinates are mapped using the offsets, scaled and
      * then cast to an integer bin.
      *
      * <pre>
      * binx = (int) (x * sx)
      * </pre>
      *
      * @param v the sample (2D coordinate xy)
      * @param observed the observations
      * @param bins the numbers of bins in the x dimension
      * @param sx the scale to convert the x coordinate to the x bin
      * @param sy the scale to convert the y coordinate to the y bin
      * @param triangle the triangle the sample should be within
      */
     private static void addObservation(double[] v, long[] observed,
             int bins, double sx, double sy, Triangle triangle) {
         final double x = v[0];
         final double y = v[1];
         // Test the point is triangle the triangle
         Assertions.assertTrue(triangle.contains(x, y));
         // Add to the correct bin after using the offset
         final int binx = (int) (x * sx);
         final int biny = (int) (y * sy);
         observed[biny * bins + binx]++;
     }
 
     /**
      * Test the SharedStateSampler implementation for 2D.
      */
     @Test
     void testSharedStateSampler2D() {
         testSharedStateSampler(2);
     }
 
     /**
      * Test the SharedStateSampler implementation for 3D.
      */
     @Test
     void testSharedStateSampler3D() {
         testSharedStateSampler(3);
     }
 
     /**
      * Test the SharedStateSampler implementation for 4D.
      */
     @Test
     void testSharedStateSampler4D() {
         testSharedStateSampler(4);
     }
 
     /**
      * Test the SharedStateSampler implementation for the given dimension.
      */
     private static void testSharedStateSampler(int dimension) {
         final UniformRandomProvider rng1 = RandomAssert.seededRNG();
         final UniformRandomProvider rng2 = RandomAssert.seededRNG();
         final double[] c1 = createCoordinate(1, dimension);
         final double[] c2 = createCoordinate(2, dimension);
         final double[] c3 = createCoordinate(-3, dimension);
         final TriangleSampler sampler1 = TriangleSampler.of(rng1, c1, c2, c3);
         final TriangleSampler sampler2 = sampler1.withUniformRandomProvider(rng2);
         RandomAssert.assertProduceSameSequence(sampler1, sampler2);
     }
 
     /**
      * Test the input vectors are copied and not used by reference for 2D.
      */
     @Test
     void testChangedInputCoordinates2D() {
         testChangedInputCoordinates(2);
     }
 
     /**
      * Test the input vectors are copied and not used by reference for 3D.
      */
     @Test
     void testChangedInputCoordinates3D() {
         testChangedInputCoordinates(3);
     }
 
     /**
      * Test the input vectors are copied and not used by reference for 4D.
      */
     @Test
     void testChangedInputCoordinates4D() {
         testChangedInputCoordinates(4);
     }
 
     /**
      * Test the input vectors are copied and not used by reference for the given
      * dimension.
      *
      * @param dimension the dimension
      */
     private static void testChangedInputCoordinates(int dimension) {
         final UniformRandomProvider rng1 = RandomAssert.seededRNG();
         final UniformRandomProvider rng2 = RandomAssert.seededRNG();
         final double[] c1 = createCoordinate(1, dimension);
         final double[] c2 = createCoordinate(2, dimension);
         final double[] c3 = createCoordinate(-3, dimension);
         final TriangleSampler sampler1 = TriangleSampler.of(rng1, c1, c2, c3);
         // Check the input vectors are copied and not used by reference.
         // Change them in place and create a new sampler. It should have different output
         // translated by the offset.
         final double offset = 10;
         for (int i = 0; i < dimension; i++) {
             c1[i] += offset;
             c2[i] += offset;
             c3[i] += offset;
         }
         final TriangleSampler sampler2 = TriangleSampler.of(rng2, c1, c2, c3);
         for (int n = 0; n < 3; n++) {
             final double[] s1 = sampler1.sample();
             final double[] s2 = sampler2.sample();
             Assertions.assertEquals(s1.length, s2.length);
             Assertions.assertFalse(Arrays.equals(s1, s2),
                     "First sampler has used the vertices by reference");
             for (int i = 0; i < dimension; i++) {
                 Assertions.assertEquals(s1[i] + offset, s2[i], 1e-10);
             }
         }
     }
 
     /**
      * Creates the coordinate of length specified by the dimension filled with
      * the given value and the dimension index: x + i.
      *
      * @param x the value for index 0
      * @param dimension the dimension
      * @return the coordinate
      */
     private static double[] createCoordinate(double x, int dimension) {
         final double[] coord = new double[dimension];
         for (int i = 0; i < dimension; i++) {
             coord[0] = x + i;
         }
         return coord;
     }
 
     /**
      * Test the transform generates 3D coordinates and can reverse them.
      */
     @Test
     void testTransform3D() {
         testTransformND(3);
     }
 
     /**
      * Test the transform generates 4D coordinates and can reverse them.
      */
     @Test
     void testTransform4D() {
         testTransformND(4);
     }
 
     /**
      * Test the transform generates ND coordinates and can reverse them.
      *
      * @param dimension the dimension
      */
     private static void testTransformND(int dimension) {
         Transform forward;
         Transform reverse;
         if (dimension == 4) {
             forward = new ForwardTransform(F4);
             reverse = new ReverseTransform(R4);
         } else if (dimension == 3) {
             forward = new ForwardTransform(F3);
             reverse = new ReverseTransform(R3);
         } else if (dimension == 2) {
             forward = reverse = new Transform2Dto2D();
         } else {
             throw new AssertionError("Unsupported dimension: " + dimension);
         }
 
         final UniformRandomProvider rng = RandomAssert.seededRNG();
         double sum = 0;
         for (int n = 0; n < 10; n++) {
             final double[] a = new double[] {rng.nextDouble(), rng.nextDouble()};
             final double[] b = forward.apply(a);
             Assertions.assertEquals(dimension, b.length);
             for (int i = 2; i < dimension; i++) {
                 sum += Math.abs(b[i]);
             }
             final double[] c = reverse.apply(b);
             Assertions.assertArrayEquals(a, c, 1e-10);
         }
         // Check that higher dimension coordinates are generated.
         // Note: Initially higher dimensions are zero.
         final double actual = sum;
         final double nonZeroThreshold = 0.01;
         Assertions.assertTrue(actual > nonZeroThreshold,
             () -> "No non-zero higher dimension coordinates: " + actual + " <= " + nonZeroThreshold);
     }
 
     /**
      * Test 3D rotations forward and reverse.
      */
     @Test
     void testRotations3D() {
         final double[] x = {1, 0.5, 0};
         final double[] y = multiply(F3, x);
         Assertions.assertArrayEquals(new double[] {0.465475314831549, 1.004183876910958, -0.157947689551155}, y, 1e-10);
         Assertions.assertEquals(length(x), length(y), 1e-10);
         final double[] x2 = multiply(R3, y);
         Assertions.assertArrayEquals(x, x2, 1e-10);
     }
 
     /**
      * Test 4D rotations forward and reverse.
      */
     @Test
     void testRotations4D() {
         final double[] x = {1, 0.5, 0, 0};
         final double[] y = multiply(F4, x);
         Assertions.assertArrayEquals(
                 new double[] {0.676776695296637, 0.780330085889911, 0.323223304703363, -0.280330085889911}, y, 1e-10);
         Assertions.assertEquals(length(x), length(y), 1e-10);
         final double[] x2 = multiply(R4, y);
         Assertions.assertArrayEquals(x, x2, 1e-10);
     }
 
     /**
      * Matrix multiplication. It is assumed the matrix is square and matches (or exceeds)
      * the vector length.
      *
      * @param m the matrix
      * @param v the vector
      * @return the result
      */
     private static double[] multiply(double[][] matrix, double[] v) {
         final int n = matrix.length;
         final int m = v.length;
         final double[] r = new double[n];
         for (int i = 0; i < n; i++) {
             double sum = 0;
             for (int j = 0; j < m; j++) {
                 sum += matrix[i][j] * v[j];
             }
             r[i] = sum;
         }
         return r;
     }
 
     /**
      * @return the length (L2-norm) of given vector.
      */
     private static double length(double[] vector) {
         double total = 0;
         for (final double d : vector) {
             total += d * d;
         }
         return Math.sqrt(total);
     }
 
     /**
      * Test the triangle contains predicate.
      */
     @Test
     void testTriangleContains() {
         final Triangle triangle = new Triangle(1, 2, 3, 1, 0.5, 6);
         // Vertices
         Assertions.assertTrue(triangle.contains(1, 2));
         Assertions.assertTrue(triangle.contains(3, 1));
         Assertions.assertTrue(triangle.contains(0.5, 6));
         // Edge
         Assertions.assertTrue(triangle.contains(0.75, 4));
         // Inside
         Assertions.assertTrue(triangle.contains(1.5, 3));
         // Outside
         Assertions.assertFalse(triangle.contains(0, 20));
         Assertions.assertFalse(triangle.contains(-20, 0));
         Assertions.assertFalse(triangle.contains(6, 6));
         // Just outside
         Assertions.assertFalse(triangle.contains(0.75, 4 - 1e-10));
 
         // Note:
         // Touching triangles can both have the point triangle.
         // This predicate is not suitable for assigning points uniquely to
         // non-overlapping triangles that share an edge.
         final Triangle triangle2 = new Triangle(1, 2, 3, 1, 0, -2);
         Assertions.assertTrue(triangle.contains(2, 1.5));
         Assertions.assertTrue(triangle2.contains(2, 1.5));
     }
 
     /**
      * Define a transform on coordinates.
      */
     private interface Transform {
         /**
          * Apply the transform.
          *
          * @param coord the coordinates
          * @return the new coordinates
          */
         double[] apply(double[] coord);
     }
 
     /**
      * Transform coordinates from 2D to a higher dimension using the rotation matrix.
      */
     private static class ForwardTransform implements Transform {
         private final double[][] r;
 
         /**
          * @param r the rotation matrix
          */
         ForwardTransform(double[][] r) {
             this.r = r;
         }
 
         @Override
         public double[] apply(double[] coord) {
             return multiply(r, coord);
         }
     }
 
     /**
      * Transform coordinates from a higher dimension to 2D using the rotation matrix.
      * The result should be in the 2D plane (i.e. higher dimensions of the transformed vector
      * are asserted to be zero).
      */
     private static class ReverseTransform implements Transform {
         private final double[][] r;
         private final int n;
 
         /**
          * @param r the rotation matrix
          */
         ReverseTransform(double[][] r) {
             this.r = r;
             n = r.length;
         }
 
         @Override
         public double[] apply(double[] coord) {
             Assertions.assertEquals(n, coord.length);
             final double[] x = multiply(r, coord);
             // This should reverse the 2D transform and return to the XY plane.
             for (int i = 2; i < x.length; i++) {
                 Assertions.assertEquals(0.0, x[i], 1e-14);
             }
             return new double[] {x[0], x[1]};
         }
     }
 
     /**
      * No-operation transform on 2D input. Asserts the input coordinates are length 2.
      */
     private static class Transform2Dto2D implements Transform {
         @Override
         public double[] apply(double[] coord) {
             Assertions.assertEquals(2, coord.length);
             return coord;
         }
     }
 
     /**
      * Class to test if a point is inside the triangle.
      *
      * <p>This function has been adapted from a StackOverflow post by Cédric Dufour. It converts the
      * point to unscaled barycentric coordinates (s,t) and tests they are within the triangle.
      * (Scaling would be done by dividing by twice the area of the triangle.)
      *
      * <h2>Warning</h2>
      *
      * <p>This assigns points on the edges as inside irrespective of edge orientation. Thus
      * back-to-back touching triangles can both have the point inside them. A predicate for geometry
      * applications where the point must be within a unique non-overlapping triangle should use
      * a different solution, e.g. assigning new points to the result of a triangulation.
      * For testing sampling within the triangle this predicate is acceptable.
      *
      * @see <a
      * href="https://stackoverflow.com/questions/2049582/how-to-determine-if-a-point-is-in-a-2d-triangle?lq=1">
      * Point in a triangle</a>
      * @see <a href="https://stackoverflow.com/a/34093754">Point inside triangle by Cédric Dufour</a>
      */
     private static class Triangle {
         private final double p2x;
         private final double p2y;
         private final double dX21;
         private final double dY12;
         private final double dY20;
         private final double dX02;
         private final double d;
 
         /**
          * Create an instance.
          *
          * @param p0 triangle vertex 0
          * @param p1 triangle vertex 1
          * @param p2 triangle vertex 2
          */
         Triangle(double[] p0, double[] p1, double[] p2) {
             this(p0[0], p0[1], p1[0], p1[1], p2[0], p2[1]);
         }
         /**
          * Create an instance.
          *
          * @param p0x triangle vertex 0 x coordinate
          * @param p0y triangle vertex 0 y coordinate
          * @param p1x triangle vertex 1 x coordinate
          * @param p1y triangle vertex 1 y coordinate
          * @param p2x triangle vertex 2 x coordinate
          * @param p2y triangle vertex 2 y coordinate
          */
         Triangle(double p0x, double p0y, double p1x, double p1y, double p2x, double p2y) {
             this.p2x = p2x;
             this.p2y = p2y;
             // Precompute factors
             dX21 = p2x - p1x;
             dY12 = p1y - p2y;
             dY20 = p2y - p0y;
             dX02 = p0x - p2x;
             // d = twice the signed area of the triangle
             d = dY12 * (p0x - p2x) + dX21 * (p0y - p2y);
         }
 
         /**
          * Check whether or not the triangle contains the given point.
          *
          * @param px the point x coordinate
          * @param py the point y coordinate
          * @return true if inside the triangle
          */
         boolean contains(double px, double py) {
             // Barycentric coordinates:
             // p = p0 + (p1 - p0) * s + (p2 - p0) * t
             // The point p is inside the triangle if 0 <= s <= 1 and 0 <= t <= 1 and s + t <= 1
             //
             // The following solves s and t.
             // Some factors are precomputed.
             final double dX = px - p2x;
             final double dY = py - p2y;
             final double s = dY12 * dX + dX21 * dY;
             final double t = dY20 * dX + dX02 * dY;
             if (d < 0) {
                 return s <= 0 && t <= 0 && s + t >= d;
             }
             return s >= 0 && t >= 0 && s + t <= d;
         }
     }
 }