Plot intersection of two filled polygons in R

Question

Let x1 and y1 be vectors of equal length, defining the coordinates of vertices of polygon1.
Let x2 and y2 be vectors of equal length, defining the coordinates of vertices of polygon2.

Polygon1, for example, can be drawn by polygon(x1,y1,border=NA,col=rgb(0,0,0))

What is the easiest way (preferably in base R i.e. without any packages) to fill only the area that belongs to both polygon1 and polygon2? In other words, what is the easiest way to draw (filled) polygon1 clipped by polygon 2?

Background:

I'm using this to shade a contoured (multi-colored) area under a standard plot. Each inter-contour region needs to be intersected with the area under the plot.

Some constraints on the polygons in my case:

In case it matters, polygon1 is in my case definable as the region between the x-axis and some y=f1(x), while polygon 2 is definable as the region enclosed between y=a*f2(x) and y=b*f2(x) where a>b.

Simplistic example data:

x1 <- 0:6
y1 <- c(0,1,2,1,0,-1,0)
x2 <- c(x1,rev(x1))
y2 <- c(x1*rev(x1)/5,x1*rev(x1)/10)

Answer 1

This draws the requisite polygon by filling everything outside the two polygons with white. Obviously it would have to be drawn before everything else in the plot, and cannot give a transparent background. If one needs to draw more than one such polygon in the same plot, this won't work.

xlim <- c(min(x1,x2),max(x1,x2))
ylim <- c(min(y1,y2),max(y1,y2))
xlim <- xlim+c(-1,1)*(xlim[2]-xlim[1])*0.04 # extend by 4% on each side
ylim <- ylim+c(-1,1)*(ylim[2]-ylim[1])*0.04 # extend by 4% on each side
n1 <- length(x1)
n2 <- length(x2)
xS <- c(xlim[1],xlim,rev(xlim),xlim[1])
yS <- c(0,rep(ylim,each=2),0)
plot(NA,xlim=xlim,ylim=ylim,xaxs='i',yaxs='i')
polygon(xS,yS,border=NA,col=rgb(0,0,0))
polygon(c(x1,xS),c(y1,yS),border=NA,col=rgb(1,1,1))
polygon(c(x2,xS),c(y2,yS),border=NA,col=rgb(1,1,1))
box()
polygon(x1,y1)
polygon(x2,y2)

So for me this answer doesn't cut it. But perhaps it could plant some seed of an idea into someone else on how to do this.

Answer 2

Since you don't give a well defined function in you example, it is impossible to define intersection region mathematically. You shoudl define the region manulaly one by one.

Another option, you can play with alpha blending to differentiate between regions.

xA <- 0:6
yA <- c(0,1,2,1,0,-1,0)
xB <- c(x1,rev(x1))
yB <- c(x1*rev(x1)/5,x1*rev(x1)/10)
plot(NA,xlim=extendrange(c(xA,xB)),
            ylim=extendrange(c(yA,yB)),xaxs='i',yaxs='i',frame.plot=F)
polygon(x = c(xA,tail(xA,1)), y = c(yA,0),col=rgb(1,0,1,0.5),border="white")
polygon(x = c(xB,tail(xB,1)), y = c(yB,0),col=rgb(1,0,1,0.5),border='white')