Affine transform with bounding box

affineBoundingBox.java

	/**
	 * The function performs an arbitrary affine transform for an arbitrary
	 * number of dimensions. The resulting view is automatically bounded
	 * 
	 * @param image
	 *            input image
	 * @param affine
	 *            defines the transformation
	 * @param interpolator
	 *            defines the type of interpolation
	 * 
	 * @return the transformed data defined on a View with bounded source
	 */
	private <K extends RealType<K>> IntervalView<K> affineBoundingBox(final RandomAccessibleInterval<K> image,
			final AffineGet affine, InterpolatorFactory<K, RandomAccessible<K>> interpolator) {

		// extend with zeros
		final RealRandomAccessible<K> field = Views.interpolate(Views.extendZero(image), interpolator);

		// transform the data
		final AffineRandomAccessible<K, AffineGet> sheared = RealViews.affine(field, affine);

		int numDims = sheared.numDimensions();

		// to hold the bounds of the transformed data
		long[] min = new long[numDims];
		long[] max = new long[numDims];
		Arrays.fill(min, Long.MAX_VALUE);
		Arrays.fill(max, Long.MIN_VALUE);

		// to hold the corners of the dataset before and after transformation
		final double[] source = new double[numDims];
		double[] target = new double[numDims];

		// dataset dimensions
		final long[] dims = new long[numDims];
		image.dimensions(dims);

		// for calculating the corner coordinates
		final long[] normDims = new long[numDims];
		Arrays.fill(normDims, 2);

		// to hold unnormalised coordinates, eg [1, 0, 0] or [1, 1, 1]
		final long[] normPos = new long[numDims];

		// loop over corners, there are 2^(numDimensions)
		for (int i = 0; i < Math.pow(2, numDims); i++) {
			// calculate unnormalised coordinates
			IntervalIndexer.indexToPosition(i, normDims, normPos);
			// normalise coordinate
			for (int d = 0; d < numDims; d++)
				source[d] = dims[d] * normPos[d];
			// apply affine transform to corner
			affine.apply(source, target);

			// if corner is extremal value in any dimension update bounds
			for (int d = 0; d < min.length; d++) {
				if (target[d] < min[d])
					min[d] = (long) Math.floor(target[d]);

				if (target[d] > max[d])
					max[d] = (long) Math.ceil(target[d]);

			}
		}
		// create interval from extremal values
		FinalInterval bounds = new FinalInterval(min, max);

		// return View with bounded source
		return Views.interval(sheared, bounds);

	}