Coronary artery disease (CAD) has one of the highest mortality rates in humans worldwide. Single-photon emission computed tomography (SPECT) myocardial perfusion imaging (MPI) provides clinicians with myocardial metabolic information non-invasively. However, there are some limitations to interpreting SPECT images performed by physicians or automatic quantitative approaches. Radiomics analyzes images objectively by extracting quantitative features and can potentially reveal biological characteristics that the human eye cannot detect. However, the reproducibility and repeatability of some radiomic features can be highly susceptible to segmentation and imaging conditions. We aimed to assess the reproducibility of radiomic features extracted from uncorrected MPI-SPECT images reconstructed with 15 different settings before and after ComBat harmonization, along with evaluating the effectiveness of ComBat in realigning feature distributions. A total of 200 patients (50% normal and 50% abnormal) including rest and stress (without attenuation and scatter corrections) MPI-SPECT images were included. Images were reconstructed using 15 combinations of filter cut-off frequencies, filter orders, filter types, reconstruction algorithms, number of iterations and subsets resulting in 6000 images. Image segmentation was performed on the left ventricle in the first reconstruction for each patient and applied to 14 others. A total of 93 radiomic features were extracted from the segmented area, and ComBat was used to harmonize them. The intraclass correlation coefficient (ICC) and overall concordance correlation coefficient (OCCC) tests were performed before and after ComBat to examine the impact of each parameter on feature robustness and to assess harmonization efficiency. The ANOVA and the Kruskal-Wallis tests were performed to evaluate the effectiveness of ComBat in correcting feature distributions. In addition, the Student's t-test, Wilcoxon rank-sum, and signed-rank tests were implemented to assess the significance level of the impacts made by each parameter of different batches and patient groups (normal vs. abnormal) on radiomic features. Before applying ComBat, the majority of features (ICC: 82, OCCC: 61) achieved high reproducibility (ICC/OCCC ≥ 0.900) under every batch except Reconstruction. The largest and smallest number of poor features (ICC/OCCC<0.500) were obtained by IterationSubset and Order batches, respectively. The most reliable features were from the first-order (FO) and gray-level co-occurrence matrix (GLCM) families. Following harmonization, the minimum number of robust features increased (ICC: 84, OCCC: 78). Applying ComBat showed that Order and Reconstruction were the least and the most responsive batches, respectively. The most robust families, in a descending order, were found to be FO, neighborhood gray-tone difference matrix (NGTDM), GLCM, gray-level run length matrix (GLRLM), gray-level size zone matrix (GLSZM), and gray-level dependence matrix (GLDM) under Cut-off, Filter, and Order batches. The Wilcoxon rank-sum test showed that the number of robust features significantly differed under most batches in the Normal and Abnormal groups. The majority of radiomic features show high levels of robustness across different OSEM reconstruction parameters in uncorrected MPI-SPECT. ComBat is effective in realigning feature distributions and enhancing radiomic features reproducibility.