---
title: "Transport plans and image interpolation"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Transport plans and image interpolation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", cache = FALSE)
library(T4transport)
```

Images are normalized pixel masses. Choose coordinates explicitly so that
comparisons and interpolation use the same geometry. The legacy image
interpolation wrapper uses normalized coordinates, increasing rightwards
and upwards. This differs from the default pixel-center coordinates of
`ot_image()`.

```{r images}
a <- b <- matrix(0, 8, 12)
a[2, 2] <- 2; a[3, 3] <- 1; a[6, 4] <- 3
b[2, 9] <- 1; b[5, 8] <- 3; b[7, 10] <- 2
x <- seq(0, 1, length.out = ncol(a))
y <- seq(1, 0, length.out = nrow(a))
mu <- ot_image(a, x = x, y = y)
nu <- ot_image(b, x = x, y = y)
exact <- ot_distance(mu, nu, method = "exact", return_plan = TRUE)
entropic <- ot_distance(mu, nu, method = "sinkhorn", epsilon = .03,
                        return_plan = TRUE)
ot_inspect(exact)
ot_inspect(entropic)
```

The entropy objective and the unregularized transport cost are separate
quantities. The rooted transport component of an entropic plan is not the
exact Wasserstein distance.

```{r interpolate}
times <- c(0, .25, .5, .75, 1)
path <- imageinterp(a, b, t = times)
lapply(path[c(1, 3, 5)], ot_inspect)
```

For each time, the returned matrix deposits mass at the nearest original
grid point. Its `unprojected_measure` attribute retains the actual
push-forward of the optimal plan. The latter follows the exact finite W2
geodesic, up to numerical precision. The displayed matrices generally do not.

```{r path-check}
mid <- attr(path[[3]], "unprojected_measure")
c(distance_to_midpoint = ot_distance(mu, mid)$wasserstein_distance,
  half_endpoint_distance = exact$wasserstein_distance / 2,
  projection_bound = attr(path[[3]], "diagnostics")$projection_w2_upper_bound)
stopifnot(max(abs(vapply(path, sum, numeric(1)) - 1)) < 1e-12)
```

```{r display, fig.width=7, fig.height=2.2}
old <- par(mfrow = c(1, 3), mar = c(1, 1, 2, 1))
for (i in c(1, 3, 5)) {
  image(x, rev(y), t(path[[i]][nrow(a):1, ]),
        axes = FALSE, asp = 1, main = paste("t =", times[i]))
}
par(old)
```

The separate replication study verifies all five times, independently
reconstructs the transported measures, and quantifies the display error.
The images here are deterministic synthetic examples with no data download.