The paper addresses the problem of virtual craniofacial reconstruction from a set of Computer Tomography (CT) images, with the multiple objectives of achieving accurate local matching of the opposable fracture surfaces and preservation of the global shape symmetry and the biomechanical stability of the reconstructed mandible. The first phase of the reconstruction, with the mean squared error as the performance metric, achieves the best possible local surface matching using the Iterative Closest Point (ICP) algorithm and the Data Aligned Rigidity Constrained Exhaustive Search (DARCES) algorithm each used individually and then in a synergistic combination. The second phase, which consists of an angular perturbation scheme, optimizes a composite reconstruction metric. The composite reconstruction metric is a linear combination of the mean squared error, a global shape symmetry term and the surface area which is shown to be a measure of biomechanical stability. Experimental results, including a thorough validation scheme on simulated fractures in phantoms of the craniofacial skeleton, are presented.