Novel graph theoretic enhancements to the well-known Iterative Closest Point (ICP) algorithm are proposed in the context of virtual craniofacial reconstruction. The input to the algorithm is a sequence of Computed Tomography (CT) images of a fractured human mandible. The closest set computation in the ICP algorithm is performed using the Maximum Cardinality Minimum Weight (MCMW) bipartite graph matching algorithm. Furthermore, the bounding boxes of the fracture surfaces are used to generate multiple candidate solutions based on the automorphism group of a cycle graph. The best candidate solution is selected by exploiting geometric constraints that are invariant to rigid body transformations and anatomical knowledge of the global shape of the mandible. Initialization of the ICP algorithm with the best candidate solution is found to improve surface reconstruction accuracy. Experimental results on CT scans of real patients are presented.