Several risk management and exotic option pricing models have been proposed in the literature which may price European options correctly. A prerequisite of these models is the interpolation of the market implied volatilities or the European option price function. However, the no-arbitrage principle places shape restrictions on the option price function. In this paper, an interpolation method is developed to preserve the shape of the option price function. The interpolation is optimal in terms of minimizing the distance between the implied risk-neutral density and the prior approximation function in L 2-norm, which is important when only a few observations are available. We reformulate the problem into a system of semismooth equations so that it can be solved efficiently.