In this article, we introduce and analyze a general procedure for approximating a ‘black and white’ shape and topology optimization problem with a density optimization problem, allowing for the presence of ‘grayscale’ regions. Our construction relies on a regularizing operator for smearing the characteristic functions involved in the exact optimization problem, and on an interpolation scheme, which endows the intermediate density regions with fictitious material properties. Under mild hypotheses on the smoothing operator and on the interpolation scheme, we prove that the features of the approximate density optimization problem (material properties, objective function, etc.) converge to their exact counterparts as the smoothing parameter vanishes. In particular, the gradient of the approximate objective functional with respect to the density function converges to either the shape or the topological derivative of the exact objective. These results shed new light on the connections between these two different notions of sensitivities for functions of the domain, and they give rise to different numerical algorithms which are illustrated by several experiments.