<p>A novel approach is presented to transform a set of coregionalized variables measured on quantitative scales into Gaussian variables. The innovation is twofold. On the one hand, unlike joint anamorphoses approaches, there is a one-to-one association between each original variable and a Gaussian variable, which makes the transformation applicable even in case of heterotopic sampling designs. On the other hand, the transformation of each variable is non-monotonic, which provides greater flexibility in comparison with the traditional normal scores transform. The inference of the transformation functions relies on the fitting of indicator direct and cross-covariances, for which an iterative procedure is proposed, and of the marginal distribution of each variable. The covariances of the Gaussian random fields can subsequently be fitted so as to reproduce the spatial correlation of (a transform of) the original variables. Once the model parameters are determined, conditional simulation can be performed by means of a mix of sequential Monte Carlo and classical multigaussian simulation techniques. An application case study is presented, pertaining to the evaluation of a lateritic nickel ore deposit, where the performances of our proposal in terms of prediction, uncertainty quantification, and reproduction of statistical and spatial dependencies between variables, are compared to that of the traditional multigaussian approach and the projection pursuit multivariate transform.</p>