From 01f07748645980b663e892366b131dee95c48586 Mon Sep 17 00:00:00 2001
From: "David C. Sterratt" <david.c.sterratt@ed.ac.uk>
Date: Mon, 31 Aug 2020 21:43:18 +0100
Subject: [PATCH] Improve image interpolation

---
 pkg/retistruct/DESCRIPTION              |  2 +-
 pkg/retistruct/NEWS                     |  2 ++
 pkg/retistruct/R/misc.R                 | 34 +++++++++++++++----------
 pkg/retistruct/man/interpolate.image.Rd |  4 ++-
 4 files changed, 26 insertions(+), 16 deletions(-)

diff --git a/pkg/retistruct/DESCRIPTION b/pkg/retistruct/DESCRIPTION
index 46bf4a67..5806231d 100644
--- a/pkg/retistruct/DESCRIPTION
+++ b/pkg/retistruct/DESCRIPTION
@@ -12,7 +12,7 @@ Description: Reconstructs retinae by morphing a flat surface with fullcuts (a
 Version: 0.7.2
 URL: http://davidcsterratt.github.io/retistruct/
 BugReports: https://github.com/davidcsterratt/retistruct/issues
-Date: 2020-08-30
+Date: 2020-08-31
 Depends:
     R (>= 3.5.0)
 Imports:
diff --git a/pkg/retistruct/NEWS b/pkg/retistruct/NEWS
index a8e4f721..b834618c 100644
--- a/pkg/retistruct/NEWS
+++ b/pkg/retistruct/NEWS
@@ -4,6 +4,8 @@ BUG FIXES
 
 * Fix regression in v0.7.1 that caused some mulit-fragment outlines to fail
 
+* Improve image interpolation
+
 CHANGES IN VERSION 0.7.1 - Released 2020/08/28
 
 NEW FEATURE
diff --git a/pkg/retistruct/R/misc.R b/pkg/retistruct/R/misc.R
index 15be04f4..958fa311 100644
--- a/pkg/retistruct/R/misc.R
+++ b/pkg/retistruct/R/misc.R
@@ -72,9 +72,10 @@ parse.dependencies <- function(deps) {
 ##' @param invert.y If \code{FALSE} (the default), the y coordinate is
 ##'   zero at the top of the image. \code{TRUE} the zero y coordinate
 ##'   is at the bottom.
+##' @param wmax maximum window size for interpolation of NA values
 ##' @return Vector of N interpolated values
 ##' @author David Sterratt
-interpolate.image <- function(im, P, invert.y=FALSE) {
+interpolate.image <- function(im, P, invert.y=FALSE, wmax=100) {
   N <- ncol(im)
   M <- nrow(im)
   x <- P[,1]
@@ -112,20 +113,25 @@ interpolate.image <- function(im, P, invert.y=FALSE) {
       t(im[c(i1,i2),c(j1,j2)]) %*%
       rbind(i1 + 0.5 - y, y - i1 + 0.5)
   }, x, y, i1, i2, j1, j2)
+
+  ## Some pixels may be in regions with NA. The stratagy here is to
+  ## extrapolate using linear regression in a window around the pixel
   for (k in which(is.na(z))) {
-    ## print(k)
-    is <- pmax(i1[k] - 10, 1):pmin(i2[k] + 10, M)
-    js <- pmax(j1[k] - 10, 1):pmin(j2[k] + 10, M)
-    dat <- data.frame(i=as.vector(outer(js*0, is, FUN="+")),
-                      j=as.vector(outer(js, is*0, FUN="+")))
-    dat$z <- mapply(function(i, j) { im[i, j] }, dat$i, dat$j)
-    dat$x  <- dat$j - 0.5
-    dat$y  <- dat$i - 0.5
-    ## print(dat)
-    ## print(dat[!is.na(dat$z),])
-    ## print(summary(lm(z ~ x + y, dat[!is.na(dat$z),c("x", "y", "z")])))
-    z[k] <- stats::predict(stats::lm(z ~ x + y, dat[!is.na(dat$z),c("x", "y", "z")]), data.frame(x=x[k], y=y[k]))
-    ## print(paste(x[k], y[k], z[k]))
+    ## Increase the window size until stats::lm() doesn't throw an error
+    w <- 1
+    while (is.na(z[k])) {
+      is <- pmax(i1[k] - w, 1):pmin(i2[k] + w, M)
+      js <- pmax(j1[k] - w, 1):pmin(j2[k] + w, M)
+      dat <- data.frame(i=as.vector(outer(js*0, is, FUN="+")),
+                        j=as.vector(outer(js, is*0, FUN="+")))
+      dat$z <- mapply(function(i, j) { im[i, j] }, dat$i, dat$j)
+      dat$x  <- dat$j - 0.5
+      dat$y  <- dat$i - 0.5
+      tryCatch(z[k] <- stats::predict(stats::lm(z ~ x + y, dat[!is.na(dat$z),c("x", "y", "z")]), data.frame(x=x[k], y=y[k])),
+               error = function(e) {}, warning = function(w) {})
+      if (w == wmax) stop(paste0("NA pixel is over ", wmax, " pixels from nearest non-NA pixel"))
+      w <- w  + 1
+    }
   }
   return(z)
 }
diff --git a/pkg/retistruct/man/interpolate.image.Rd b/pkg/retistruct/man/interpolate.image.Rd
index a9c2f1dd..a6041c98 100644
--- a/pkg/retistruct/man/interpolate.image.Rd
+++ b/pkg/retistruct/man/interpolate.image.Rd
@@ -4,7 +4,7 @@
 \alias{interpolate.image}
 \title{Interpolate values in image}
 \usage{
-interpolate.image(im, P, invert.y = FALSE)
+interpolate.image(im, P, invert.y = FALSE, wmax = 100)
 }
 \arguments{
 \item{im}{image to interpolate}
@@ -15,6 +15,8 @@ is in range \code{[0, ncol(im)]} and y is in range \code{[0, nrow(im)]}}
 \item{invert.y}{If \code{FALSE} (the default), the y coordinate is
 zero at the top of the image. \code{TRUE} the zero y coordinate
 is at the bottom.}
+
+\item{wmax}{maximum window size for interpolation of NA values}
 }
 \value{
 Vector of N interpolated values