Improve image interpolation

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:

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

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)
 }