In-silico modeling of atrial fibrillation (AF) requires anatomically realistic, population-representative shape models that integrate with electrophysiology (EP) modeling frameworks. However, existing statistical shape models (SSMs) often rely on small cohorts or provide limited anatomical coverage. We constructed left atrial (LA), right atrial (RA), and bi-atrial SSMs from late gadolinium-enhanced MRI scans of 74 AF patients. The pipeline involved manual segmentation, rigid alignment, registration, and principal component analysis using the open-source library Scalismo (Scalable Image Analysis and Shape Modeling). From each SSM, 1000 synthetic meshes were generated by sampling principal components (PCs) within ±3 standard deviations. Atrial fiber orientations from a bilayer atlas were mapped onto the synthetic geometries using the atrialmtk framework. This pipeline generated 1000 personalized bi-atrial bilayer meshes for use in EP simulations to assess arrhythmia dynamics and guide patient-specific treatment strategies. The bi-atrial SSM required 29 PCs to capture 95% of shape variance, compared to 21 for the LA-only and 24 for the RA-only SSM. The leading modes represented variation in atrial size, RA sphericity, and pulmonary vein configuration. The synthetic meshes had physiologically realistic volumes, surface areas, and low specificity RMSE (5.82±0.53 mm). Finite element simulations were performed and showed AF was inducible in the synthetic bi-atrial models. The bi-atrial synthetic models captured inter-chamber morphological coupling and fiber continuity better than independently generated LA and RA models. This pipeline provides a scalable approach for generating personalized bi-atrial bilayer meshes for EP studies.