Reliability-based design optimization (RBDO) algorithms have been developed to solve design optimization problems with existence of uncertainties. Traditionally, the original random design space is transformed to the standard normal design space, where the reliability index can be measured in a standardized unit. In the standard normal design space, the modified reliability index approach (MRIA) measured the minimum distance from the design point to the failure region to represent the reliability index; on the other hand, the performance measure approach (PMA) performed inverse reliability analysis to evaluate the target function performance in a distance of reliability index away from the design point. MRIA was able to provide stable and accurate reliability analysis while PMA showed greater efficiency and was widely used in various engineering applications. However, the existing methods cannot properly perform reliability analysis in the standard normal design space if the transformation to the standard normal space does not exist or is difficult to determine. To this end, a new algorithm, ensemble of Gaussian reliability analyses (EoGRA), was developed to estimate the failure probability using Gaussian-based kernel density estimation (KDE) in the original design space. The probabilistic constraints were formulated based on each kernel reliability analysis for the optimization processes. This paper proposed an efficient way to estimate the constraint gradient and linearly approximate the probabilistic constraints with fewer function evaluations (FEs). Some numerical examples with various random distributions are studied to investigate the numerical performances of the proposed method. The results showed that EoGRA is capable of finding correct solutions in some problems that cannot be solved by traditional methods. Furthermore, experiments of image processing with arbitrarily distributed photo pixels are performed. The lighting of image pixels is maximized subject to the acceptable limit. Our implementation showed that the accuracy of the estimation of normal distribution is poor while the proposed method is capable of finding the optimal solution with acceptable accuracy.