The basic problem of shape complementarity analysis appears fundamental to applications as diverse as mechanical design, assembly automation, robot motion planning, micro- and nano-fabrication, protein-ligand binding, and rational drug design. However, the current challenge lies in the lack of a general mathematical formulation that applies to objects of arbitrary shape. We propose that a measure of shape complementarity can be obtained from the extent of approximate overlap between shape skeletons. A space-continuous implicit generalization of the skeleton, called the skeletal density function (SDF) is defined over the Euclidean space that contains the individual assembly partners. The SDF shape descriptors capture the essential features that are relevant to proper contact alignment, and are considerably more robust than the conventional explicit skeletal representations. We express the shape complementarity score as a convolution of the individual SDFs. The problem then breaks down to a global optimization of the score over the configuration space of spatial relations, which can be efficiently implemented using fast Fourier transforms (FFTs) on nonequispaced samples. We demonstrate the effectiveness of the scoring approach for several examples from 2D peg-in-hole alignment to more complex 3D examples in mechanical assembly and protein docking. We show that the proposed method is reliable, inherently robust against small perturbations, and effective in steering gradient-based optimization.