A procedure to reconstruct two-phase porous media, given the porosity and the two-point correlation function of such media, is described. The random media are modelled as a discrete valued random field Z(x), where x is the position vector. Z(x) takes value 1 in regions of pores and 0 in regions of solid phase. The field Z(x) is obtained by applying a non-linear filter—Nataf’s transformation—to a correlated Gaussian random field Y(x). The two-point correlation function R YY of the Gaussian field Y is related to the two-point correlation function R ZZ of the field Z and can be calculated by expanding the bivariate Gaussian probability density in terms of Hermite polynomials. This expansion results in a series representation for R ZZ in terms of R YY . Depending on the accuracy intended, the series could be truncated, and the appropriate root of the polynomial equation thus obtained gives R YY . The correlation function of the Gaussian field is decomposed into eigenfunctions and eigenvalues required by the Karhunen-Lóeve expansion. The eigenfunctions and eigenvalues are used to generate as many samples of the Gaussian field as required, and the discrete field corresponding to each such sample can be obtained by applying the non-linear filter mentioned above. The method was tested by generating a large number of samples of two-dimensional Debye random media using different porosities and different correlation lengths, and the statistics of the ensemble was found to agree favourably with the input data. The method also has the advantage that it gives a theoretical framework for the porous media in terms of the random fields. These random fields could be used to model fluid flow through such porous media.