In this work, we propose a Bayesian nonparametric scalar-on-image regression model, where the output is a scalar and the image is treated as the input, with a thresholded Gaussian process prior on the spatially-varying coefficients in order to introduce sparsity and smoothness into the model. Specifically, we apply a class of piecewise-smooth, sparse, and continuous functions, called the relaxed-thresholded Gaussian process (RTGP), that are flexible enough to exhibit both properties of hard- and soft-thresholding Gaussian processes to tackle the common non-identifiability issue of scalar-on-image regression problems. Our main contribution is the improved scalability, allowing for larger sample sizes and bigger image dimensions, which is made possible by replacing posterior sampling with a variational approximation. Moreover, the reduction in computational cost enables the vertex-wise analysis of cortical surface data or the voxel-wise analysis of volumetric data which previously was too prohibitive to perform for large-scale Bayesian spatial models with thresholded Gaussian process priors. Lastly, we validate our results via extensive simulation studies and an application to the Adolescent Brain Cognition Development (ABCD) study, a large-scale neuroimaging application with a sample size of over 5,000 subjects.