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→), which 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 RYY of the Gaussian field Y is related to the two point correlation function RZZ of the field Z and can be calculated by expanding the bivariate Gaussian probability in terms of Hermite polynomials. 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 one and 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. Also one and two dimensional 'chess board' patterns were reconstructed to see how well the geometry is reconstructed. The one dimensional case was reconstructed very accurately, whereas the two dimensional case, though not very satisfactory, indicates that the method captures some of the essentials of the geometry. The method also has the advantage that it gives an analytical framework for the porous media in terms of the random fields. These random fields could be used for further studies related to porous media. © 2013 Elsevier Ltd.