This is potentially a very difficult problem when the borders are ragged. A brute-force search of the optimum could require computational time that is proportional to the square of the number of cells in the image (a value that often will be in the billions, trillions, or greater).
One promising approach is to relax the conditions a little bit and actually allow a few NoData values to sneak in: simply penalize their inclusion. This affords some welcome flexibility, too: a small string of NoData values penetrating from the outside (such as a skinny gap) could be included rather than cutting off most of the image.
This suggests maximizing the area of the remaining image, subject to a penalty that depends on how many NoData cells are left in the remaining image. To solve this problem it helps to interpolate positions within the image, so that it changes from a large discrete problem (where only integral row and column indexes are allowed) to a continuous problem. Bilinear interpolation should be fine.
Only four variables are needed to represent the coordinates of the lower left and upper right of the final included image. This makes the problem quite tractable.
To expedite the calculations of the numbers of included NoData cells, first compute a cumulative sum array Y by summing along columns and then accumulating those sums along rows. To obtain the number of NoData cells included in the rectangle spanning the coordinates (i,k) through (j,l), just calculate
Y[k,l] - Y[k,j] - Y[i,l] + Y[i,j]
That's much faster than re-summing them each time.
A canned optimizer (in R) had some difficulties finding solutions. I was able to tweak it by starting with a small penalty factor and gradually increasing it, restarting each new search at the value of the preceding one. Occasionally this fails and needs to be restarted with a different value, but most of the time it works beautifully.
Because of the interpolation, an extra row or column of pixels might be stripped away unnecessarily. A little post-processing could detect this circumstance and fix it up. (I haven't implemented that.)
Here is the progress of one calculation involving an 85 by 119 array with randomly generated NoData cells within 15 pixels of its border all around. The upper row plots the masked image (starting with the original at the upper left). The bottom row, from left to right, depicts the cumulative sum array Y and the successive masks that were found. The final one includes 4717 pixels, which is slightly less than the correct value of 4895 pixels (because of the aforementioned interpolation issue). This calculation took 0.3 seconds.
It scales beautifully: on a similar configuration with 100 times as many pixels the optimization step took under one second, only three times longer. Another factor of 100 (an image of over 100,000,000 pixels) brought the time up to 80 seconds. (You have to turn off the display of the intermediate images--R is particularly slow at that.)
The R code that produced this example is a bit of a hack--it's not bulletproof nor is it optimized for speed. It is offered to illustrate the details and perhaps to serve as a point of departure for anyone who would like to improve or port it.
#
# Create data.
#
border <- 15
m <- 55+2*border
n <- 89+2*border
set.seed(17)
x <- rep(0, m*n)
e <- runif(m*n) < 1/3
i <- as.vector(outer(1:m, 1:n,
function(a,b) (a < 1+border | a > m-border) |
(b < 1+border | b > n-border)))
x[i] <- e[i]
x <- matrix(x, m, n)
par(mfcol=c(2,4))
#
# Display the original.
#
n.bad <- sum(x)
image(x, col=c("Gray", "Red"), main=paste("Bad pixels =", n.bad))
#
# Compute the cumulative sums of NA data.
#
y <- apply(apply(x, 1, cumsum), 1, cumsum)
image(y, col=rainbow(max(m,n)), main="Cumulative bad pixels")
#
# Perform bilinear interpolation into an array.
#
access <- function(y, ij) {
bounds <- function(i, n) {
lower <- pmin(n-1, pmax(1, floor(i)))
upper <- lower+1
delta <- pmax(0, pmin(1, c(upper-i, i-lower)))
return (list(lower=lower, upper=upper, delta=delta))
}
b.1 <- bounds(ij[1], dim(y)[1])
b.2 <- bounds(ij[2], dim(y)[2])
values <- c(y[b.1$lower,b.2$lower], y[b.1$upper,b.2$lower],
y[b.1$lower,b.2$upper], y[b.1$upper,b.2$upper])
weights <- as.vector(outer(b.1$delta, b.2$delta, '*'))
return (sum(values*weights))
}
# z <- matrix(NA, floor(pi*m), floor(pi*n))
# for (i in 1:dim(z)[1]) for (j in 1:dim(z)[2]) z[i,j] <- access(y, c(i,j)/pi)
# image(z)
#
# Define the objective function (to be minimized).
#
f <- function(indexes, y, rho=1, rho.boundary=1) {
clamp <- function(x, d) pmin(d, pmax(1, x))
m <- dim(y)[1]; n <- dim(y)[2]
i0 <- min(indexes[1:2])
k0 <- max(indexes[1:2])
j0 <- min(indexes[3:4])
l0 <- max(indexes[3:4])
i <- clamp(i0, m)
k <- clamp(k0, m)
j <- clamp(j0, n)
l <- clamp(l0, n)
area <- abs((k-i+1) * (l-j+1))
penalty <- access(y, c(k,l)) - access(y, c(k,j)) - access(y, c(i,l)) + access(y, c(i,j))
penalty.b <- sum((c(i,j,k,l) - c(i0,j0,k0,l0))^2)
return (-area + rho*penalty + rho.boundary*penalty.b)
}
# f.r <- matrix(NA, floor(pi*m), floor(pi*n))
# for (i in 1:dim(z)[1]) for (j in 1:dim(z)[2]) f.r[i,j] <- f(c(15,i,7,j), z)
# image(f.r)
#
# Sneak up on a solution by starting with small penalties and
# increasing them gradually.
#
theta.0 <- c(1,m, 1,n)
rho <- 1/2
n.max <- 16
n.bad.old <- n.bad
repeat {
if (n.max <= 0) break
sol <- nlm(function(theta) f(theta, y, rho, 1), theta.0, steptol=0.1)
mask <- matrix(FALSE, m, n)
ij <- sol$estimate
i0 <- ceiling(min(ij[1:2]))
k0 <- floor(max(ij[1:2]))
j0 <- ceiling(min(ij[3:4]))
l0 <- floor(max(ij[3:4]))
mask[i0:k0, j0:l0] <- TRUE
x.mask <- x
x.mask[!mask] <- NA
#image(x)
n.bad <- sum(x.mask, na.rm=TRUE)
if (n.bad < n.bad.old) {
image(x.mask, col=c("Gray", "Red"),
main=paste("Bad pixels =", n.bad), sub=paste("Penalty =", round(rho, 1)))
image(mask, col=c("Orange", "White"), main="Mask")
n.bad.old <- n.bad
}
if (n.bad == 0) break
n.max <- n.max - 1
rho <- rho*1.5
theta.0 <- sol$estimate
}
#
# Compare the number of displayed pixels to the best value.
#
sum(mask)
(m-2*border)*(n-2*border)
