Tortuosity is a parameter that plays a significant role in the characterization of complex porous media systems and it has a significant impact on many engineering and environmental processes and applications. Flow in porous media, diffusion of gases in complex pore structures and membrane flux in water desalination are examples of the application of this important micro-scale parameter. In this paper, an algorithm was developed and implemented as a MATLAB code to compute tortuosity from three-dimensional images. The code reads a segmented image and finds all possible tortuous paths required to compute tortuosity. The code is user-friendly, easy to use and computationally efficient, as it requires a relatively short time to identify all possible connected paths between two boundaries of large images. The main idea of the developed algorithm is that it conducts a guided search for connected paths in the void space of the image utilizing the medial surface of the void space. Once all connected paths are identified in a specific direction, the average of all connected paths in that direction is used to compute tortuosity. Three-dimensional images of sand systems acquired using X-ray computed tomography were used to validate the algorithm. Tortuosity values were computed from three-dimensional images of nine different natural sand systems using the developed algorithm and compared with predicted values by models available in the literature. Findings indicate that the code can successfully compute tortuosity for any unconsolidated porous system irrespective of the shape (i.e., geometry) of particles.(C) 2017 Elsevier B.V. All rights reserved.