Lipid membrane packing defects are considered to be an essential parameter that regulates specific membrane binding of several peripheral proteins. In the absence of direct experimental characterization, lipid packing defects and their role in the binding of peripheral proteins are generally investigated through computational studies, which have been immensely successful in unraveling the key steps of the membrane-binding process. However, packing defects are calculated using two-dimensional (2D) projections and the crucial information on their depths is generally overlooked. Here, we present a simple yet computationally efficient algorithm, which identifies these defects in three dimensions. We validate the algorithm on a number of model membrane systems that are previously studied using 2D defect calculations and find that the defect size and the defect depth may not always be directly correlated. We employ the algorithm to understand the nature of packing defects in flat bilayer membranes exhibiting liquid-ordered (Lo), liquid-disordered (Ld), and co-existing (Lo/Ld) phases. Our results indicate the presence of shallower, smaller, and spatially localized defects in the Lo phase membranes as compared to the defects in Ld and mixed Lo/Ld phase membranes. Such analyses can elucidate the molecular-scale mechanisms that drive the preferential localization of certain proteins to either of the liquid phases or their interface. We also analyze the membrane sensing and anchoring process of a peptide in terms of the three-dimensional defects, which provides additional insights into the process with respect to depth distributions across the bilayer leaflets.