Motivated by recent advances in technology for medical imaging and high-throughput genotyping, we consider an imaging genetics approach to discover relationships between the interplay of genetic variation and environmental factors and measurements from imaging phenotypes. We propose an image-on-scalar regression method, in which the spatial heterogeneity of gene-environment interactions on imaging responses is investigated via an ultra-high-dimensional spatially varying coefficient model (SVCM). Bivariate splines on triangulations are used to represent the coefficient functions over an irregular two-dimensional (2D) domain of interest. For the proposed SVCMs, we further develop a unified approach for simultaneous sparse learning (i.e., G×E interaction identification) and model structure identification (i.e., determination of spatially varying vs. constant coefficients). Our method can identify zero, nonzero constant and spatially varying components correctly and efficiently. The estimators of constant coefficients and varying coefficient functions are consistent and asymptotically normal. The performance of the method is evaluated by Monte Carlo simulation studies and a brain mapping study based on the Alzheimer's Disease Neuroimaging Initiative (ADNI) data.