Nonparametric estimation: Interaction of two continuous variables The probability of observing a binary outcome y is a function of two variables, y = f(x1, x2).

Both x1 and x2 are continuous, and I have no prior on the functional form of f. I have shallow knowledge of kernel density estimation, but have only seen that done in the case of a univariate variable. How should I proceed here / is there some basic literature on this?
(I misunderstood the question at the beginning) 
The model is a binary non parametric regression. 
I don't know of a general model, I know of a semiparametric model by Klein and Spady 
http://www.ssc.wisc.edu/~bhansen/718/NonParametrics7.pdf
The model is as follows let $x´=[x_1´ \quad x_2´]$ 
Now $P(y=1|x)=E[y=1|x]=g(x'\beta)=f(x)$. 
The details in the link. 
Another possibility is to estimate the joint distribution (that has both discrete and continuous variables). Let g(y,x) be the joint pdf, then estimate $g(y|x)=\frac{g(y,x)}{g(x)}=f(x)$. The estimation of a joint of this type needs a product kernel that you can check 
Nonparametric estimation of regression functions with both categorical and continuous data
Journal of Econometrics, Volume 119, Issue 1, March 2004, Pages 99–130
Jeff Racine, Qi Li

For $y$ continuous the following answer works.
First of all is your data deterministic or random. 
If deterministic you may wish to do interpolation. Splines are your best option. 
If you think your dataset has noise or you have unobserved variables then you can do a non-parametric regression, I think this is what you want. 
The first step is to actually make assumptions about the structure of the randomness. 
The usual assumption is that 
$y=f(x_1,x_2)+\epsilon$ and you let $x´=[x_1´ \quad x_2´]$ be orthogonal -indepdent- of the unobserved variable $x\perp\epsilon$. Then you can replace $f$ by some high degree polynomial on $x$ and estimate this using Non linear least squares, technically this is called series or polynomial estimators (NLS is used when you know f, but you can argue that a high degree polynomial approximates the true f).
Another way is using a ¨kernel¨ approach. 
In fact  you can directly assume that f is the conditional expectation. 
$E[y|x]=f(x_1,x_2)$. In that case you can estimate it using the joint distribution recall the definition of conditional expectation 
$E[y|x]=\int yf(y|x)dy$, with $f(y|x)=\frac{f(y,x)}{f(x)}$. Replacing this elements by its kernel estimators you get the Nadaraya-Watson estimator, you can google it. 
A third option is to use local linear regression. A more advanced option is to use Sieves estimators that use function basis to estimate the model. 
If f is not the conditional expectation then you have to use non parametric regression with instrumental variables, this is more complicated there is a paper by Renault (and others) on Econometrica about this. 
In summary it depends on the randomness structure that you are assuming. 
The book to go if you want to learn non-parametric regression is Li &amp; Racine http://press.princeton.edu/titles/8355.html
There is also a nice implementation of these econometric procedures in the open source statistical language "R" under the name (package) "np".