[R] shading in image()

[Angel]

Is there a way to make a shading interpolation on an image plot?
Something similar to matlab 'shading interp', I think it is called Gouraud
shading.

What I want is to make a image plot look nicer. with image() it looks very
facetted, and I would like to make it look smoother.
I've tried with interp.surface() in fields package but it (obviously) makes
nan values at the borders and around nan.

Maybe an example would make me explain better:

data(volcano)
par(mfrow=c(1,2))
x <- 10*(1:nrow(volcano))
y <- 10*(1:ncol(volcano))
volcano[seq(1,87,by=10),seq(1,61,by=10)]<-NA;
image(x, y, volcano, col = terrain.colors(100))
## This is what I've tried with interp.surface()
library(fields)
newx<-seq(10,870,by=2);
newy<-seq(10,610,by=2);
volcano2<-interp.surface(list(x=x,y=y,z=volcano),make.surface.grid(list(newx
,newy)));
dim(volcano2)<-c(length(newx),length(newy));
image(newx,newy,volcano2, col = terrain.colors(100))

As you see NANs are lareger in the interpolated image, I've modified the
interp.surface function to take the nearest neightbour when any of the 4
bounding points is NA but the result is still not perfect, thats why I'm
looking for a function that does Gouraud shading.
Thanks for all the help,
Angel