How to mask or crop IDW results
Here I am answering my own question:

To obtain an IDW surface interpolation restricted to the area containing the original set of xyz data points:


Create a rectangular grid based on the xy point ranges using expand.grid {base} 
Create a convex hull using the function convexhull.xy{spatstat}
Cut an irregular grid using the function inout{splancs} on the
rectangular grid masked by the convex hull
Perform the interpolation using the function idw{gstat} with the new grid.
Be sure to detach spatstat, because it has a different function for
IDW with the same name


Code:

###
# IDW in R
#
# Load gstat:
library(gstat)
# 
# An example XYZ data set:
set.seed(1234)
x &lt;- rnorm(20)*1000+3000
y &lt;- rnorm(20)*1000+3000
z &lt;- rnorm(20)+10
xyz &lt;- data.frame(x,y,z)
# 
# Calculate ranges and spacing:
xr &lt;- round(range(x),1);xr
yr &lt;- round(range(y),1);yr
b &lt;- round((xr[2]-xr[1])/60,2);b
#
# Create a grid:
grd1 &lt;- expand.grid(x=seq(from=xr[1],to=xr[2],by=b),
    y=seq(from=yr[1],to=yr[2],by=b))
plot(grd1)
#
# Load spatstat and splancs
library(spatstat)
library(splancs)
#
# Extract convex hull 
w2 &lt;- convexhull.xy(x,y)
#
# detach spatstat because it has a different "idw" function:
detach("package:spatstat", unload=TRUE)
#
# Create a new cropped grid:
grd2 &lt;- grd1[inout(grd1,w2$bdry[[1]]),]
plot(grd2)
#
# Convert grids and dataset to sp objects:
coordinates(grd1) &lt;- ~x+y
gridded(grd1)   &lt;- TRUE
coordinates(grd2) &lt;- ~x+y
gridded(grd2)   &lt;- TRUE
coordinates(xyz) &lt;- ~x+y
#
# Run the interpolations:
idw1&lt;-idw(formula=z~1, locations=xyz, newdata=grd1)
idw2&lt;-idw(formula=z~1, locations=xyz, newdata=grd2)
#
# Plot the results:
image(idw1, col=terrain.colors(16))
contour(idw1, add=T)
points(x,y, col=4, pch=19)
image(idw2, col=terrain.colors(16))
contour(idw2, add=T)
points(x,y, col=4, pch=19)
#