Cellular mechanics drives epithelial morphogenesis, the process wherein cells collectively rearrange to produce tissue-scale deformations that determine organismal shape. However, quantitative understanding of tissue mechanics is impaired by the difficulty of direct measurement of stress in-vivo. This difficulty has spurred the development of image-based inference algorithms that estimate stress from snapshots of epithelial geometry. Such methods are challenged by sensitivity to measurement error and thus require accurate geometric segmentation for practical use. We overcome this difficulty by introducing a novel approach - the Variational Method of Stress Inference (VMSI) - which exploits the fundamental duality between stress and geometry at equilibrium of discrete mechanical networks that model confluent cellular layers. We approximate the apical geometry of an epithelial tissue by a 2D tiling with Circular Arc Polygons (CAP) in which arcs represent intercellular interfaces defined by the balance of local line tension and pressure differentials between adjacent cells. The mechanical equilibrium of such networks imposes extensive local constraints on CAP geometry. These constraints provide the foundation of VMSI which, starting with images of epithelial monolayers, simultaneously approximates both tissue geometry and internal forces, subject to the constraint of equilibrium. We find VMSI to be more robust than previous methods. Specifically, the VMSI performance is validated by the comparison of the predicted cellular and mesoscopic scale stress with the measured myosin II patterns during early Drosophila embryogenesis. VMSI prediction of mesoscopic stress tensor correlates at the 80% level with the measured myosin distribution and reveals that most of the myosin activity in that case is involved in a static internal force balance within the epithelial layer. In addition to insight into cell mechanics, this study provides a practical method for non-destructive estimation of stress in live epithelial tissue.