Nonparametric estimation of copula density functions using kernel estimators presents significant challenges. One issue is the potential unboundedness of certain copula density functions at the corners of the unit square. Another is the boundary bias inherent in kernel density estimation. This paper presents a kernel-based method for estimating bounded copula density functions, addressing boundary bias through the mirror-reflection technique. Optimal smoothing parameters are derived via Asymptotic Mean Integrated Squared Error (AMISE) minimization and cross-validation, with theoretical guarantees of consistency and asymptotic normality. Two kernel smoothing strategies are proposed: the rule-of-thumb approach and least squares cross-validation (LSCV). Simulation studies highlight the efficacy of the rule-of-thumb method in bandwidth selection for copulas with unbounded marginal supports. The methodology is further validated through an application to the Wisconsin Breast Cancer Diagnostic Dataset (WBCDD), where LSCV is used for bandwidth selection.