The following example assumes a stationary point process.
library(spatstat)
data(simdat)
X = simdat
plot(X)

E = envelope(X, Kest,nsim=99)
plot(E, sqrt(./pi)-r ~ r)

This example does NOT assume a stationary process.
library(spatstat)
data(lansing)
# inhomogeneous pattern of maples
X = unmark(split(lansing)$maple)
plot(X)

# Intensity function estimated by model-fitting
# Fit spatial trend: polynomial in x and y coordinates
fit = ppm(X, ~ polynom(x,y,2), Poisson())
# (a) predict intensity values at points themselves,
# obtaining a vector of lambda values
lambda = predict(fit, locations=X, type="trend")
# inhomogeneous K function
Ki = Kinhom(X, lambda)
plot(Ki)

# (b) predict intensity at all locations,
# obtaining a pixel image
lambda = predict(fit, type="trend")
plot(lambda)
plot(X,add=T)

Ki = Kinhom(X, lambda)
plot(Ki)

# How to make simulation envelopes:
# Example shows method (1b)
lambda = predict(fit, type="trend")
Ken = envelope(X, Kinhom, nsim=99,
simulate=expression(rpoispp(lambda)),
lambda=lambda, correction="trans")
plot(Ken, sqrt(.) - r ~ r)
