Pore network models, in which the pore space is represented by a 3D network of interconnected pores and throats, are used extensively to compute important macroscopic transport properties including capillary pressure, relative permeability and residual saturation [2,6,12]. The predictive value of network models depends on the accuracy with which the network captures the complex geometric and topological properties of real porous rocks. A practical approach for a wide range of rock types is to extract networks and network properties directly from high-resolution 3D images of the pore space [2,7,8]. To ensure that generated networks are accurate representations of the imaged rock one must overcome problems of sensitivity to image noise and the lack of a robust procedure for merging adjacent nodes to form pores. One must also be able to generate networks on 3D volumes that are sufficiently large to be representative . We present an evolutionary approach for network extraction that uses the medial axis transform together with a number of morphological measures to select tessellation boundaries and applies a new node merging algorithm. The algorithms are fully parallel, allowing very large networks containing up to a million nodes to be generated. The power and flexibility of the network extraction procedure is illustrated by examining micro-CT images for a number of sandstone and carbonate samples at image sizes of up to 2000 3 voxels and resolutions down to 2 microns. The variability in network structure obtained across the range of samples imaged highlights the need to generate realistic pore network structures when attempting to perform predictive two phase flow modeling.