	@Test
	public void testSamplePoints() {

		// SETUP
		final double a = 2.0;
		final double b = 3.0;
		final double c = 4.0;
		final Ellipsoid ellipsoid = new Ellipsoid(a, b, c);
		final Vector3dc centroid = new Vector3d(4, 5, 6);
		ellipsoid.setCentroid(centroid);
		final double sinAlpha = Math.sin(Math.PI / 4.0);
		final double cosAlpha = Math.cos(Math.PI / 4.0);
		// Orientation of an ellipsoid that's been rotated 45 deg around the z-axis
		// @formatter:off
		final Matrix3dc orientation = new Matrix3d(
				cosAlpha, -sinAlpha, 0,
				sinAlpha, cosAlpha, 0,
				0, 0, 1
		);
		// @formatter:on
		final AxisAngle4d invRotation = new AxisAngle4d();
		invRotation.set(orientation);
		invRotation.angle = -invRotation.angle;
		final Quaterniondc q = new Quaterniond(new AxisAngle4d(invRotation.angle, invRotation.x, invRotation.y, invRotation.z));
		final BinaryHybridCFI1<Vector3d, Quaterniondc, Vector3d> inverseRotation =
			Hybrids.binaryCFI1(IMAGE_J.op(), Rotate3d.class, Vector3d.class,
					new Vector3d(), q);
		ellipsoid.setOrientation(orientation);
		final long n = 10;
		final Function<Vector3d, Double> ellipsoidEq = (Vector3d p) -> {
			final BiFunction<Double, Double, Double> term = (x, r) -> (x * x) / (r *
				r);
			return term.apply(p.x, a) + term.apply(p.y, b) + term.apply(p.z, c);
		};
		ellipsoid.initSampling(IMAGE_J.op());

		// EXECUTE
		final List<Vector3d> points = ellipsoid.samplePoints(n);

		// VERIFY
		assertNotNull(points);
		assertEquals(n, points.size());
		// Reverse the translation and rotation so that we can assert the easier
		// equation of an axis-aligned ellipsoid at centroid
		points.forEach(p -> p.sub(centroid));
		points.forEach(inverseRotation::mutate);
		points.forEach(p -> assertEquals(
			"Point doesn't solve the ellipsoid equation", 1.0, ellipsoidEq.apply(p),
			1e-5));
	}