diff --git a/data/.gitattributes b/data/.gitattributes
index faa46af5..792ce382 100644
--- a/data/.gitattributes
+++ b/data/.gitattributes
@@ -1,2 +1,3 @@
 *.vtk filter=lfs diff=lfs merge=lfs -text
 *.xyz filter=lfs diff=lfs merge=lfs -text
+*.ply filter=lfs diff=lfs merge=lfs -text
diff --git a/data/cube.ply b/data/cube.ply
new file mode 100644
index 00000000..a4ec5530
--- /dev/null
+++ b/data/cube.ply
@@ -0,0 +1,51 @@
+ply
+format ascii 1.0
+comment Created by Blender 2.82 (sub 7) - www.blender.org, source file: ''
+element vertex 24
+property float x
+property float y
+property float z
+property float nx
+property float ny
+property float nz
+property float s
+property float t
+element face 12
+property list uchar uint vertex_indices
+end_header
+-1.000000 1.000000 1.000000 0.000000 0.000000 1.000000 0.875000 0.500000
+1.000000 -1.000000 1.000000 0.000000 0.000000 1.000000 0.625000 0.750000
+1.000000 1.000000 1.000000 0.000000 0.000000 1.000000 0.625000 0.500000
+1.000000 -1.000000 1.000000 0.000000 -1.000000 0.000000 0.625000 0.750000
+-1.000000 -1.000000 -1.000000 0.000000 -1.000000 0.000000 0.375000 1.000000
+1.000000 -1.000000 -1.000000 0.000000 -1.000000 0.000000 0.375000 0.750000
+-1.000000 -1.000000 1.000000 -1.000000 0.000000 0.000000 0.625000 0.000000
+-1.000000 1.000000 -1.000000 -1.000000 0.000000 0.000000 0.375000 0.250000
+-1.000000 -1.000000 -1.000000 -1.000000 0.000000 0.000000 0.375000 0.000000
+1.000000 1.000000 -1.000000 0.000000 0.000000 -1.000000 0.375000 0.500000
+-1.000000 -1.000000 -1.000000 0.000000 0.000000 -1.000000 0.125000 0.750000
+-1.000000 1.000000 -1.000000 0.000000 0.000000 -1.000000 0.125000 0.500000
+1.000000 1.000000 1.000000 1.000000 0.000000 -0.000000 0.625000 0.500000
+1.000000 -1.000000 -1.000000 1.000000 0.000000 -0.000000 0.375000 0.750000
+1.000000 1.000000 -1.000000 1.000000 0.000000 -0.000000 0.375000 0.500000
+-1.000000 1.000000 1.000000 0.000000 1.000000 -0.000000 0.625000 0.250000
+1.000000 1.000000 -1.000000 0.000000 1.000000 -0.000000 0.375000 0.500000
+-1.000000 1.000000 -1.000000 0.000000 1.000000 -0.000000 0.375000 0.250000
+-1.000000 -1.000000 1.000000 0.000000 -0.000000 1.000000 0.875000 0.750000
+-1.000000 -1.000000 1.000000 0.000000 -1.000000 0.000000 0.625000 1.000000
+-1.000000 1.000000 1.000000 -1.000000 0.000000 0.000000 0.625000 0.250000
+1.000000 -1.000000 -1.000000 0.000000 0.000000 -1.000000 0.375000 0.750000
+1.000000 -1.000000 1.000000 1.000000 0.000000 0.000000 0.625000 0.750000
+1.000000 1.000000 1.000000 0.000000 1.000000 -0.000000 0.625000 0.500000
+3 0 1 2
+3 3 4 5
+3 6 7 8
+3 9 10 11
+3 12 13 14
+3 15 16 17
+3 0 18 1
+3 3 19 4
+3 6 20 7
+3 9 21 10
+3 12 22 13
+3 15 23 16
diff --git a/splashsurf/src/io.rs b/splashsurf/src/io.rs
index 4e379f18..a3a808a0 100644
--- a/splashsurf/src/io.rs
+++ b/splashsurf/src/io.rs
@@ -143,6 +143,7 @@ pub fn read_surface_mesh<R: Real, P: AsRef<Path>>(
 
         match extension.to_lowercase().as_str() {
             "vtk" => vtk_format::surface_mesh_from_vtk(&input_file)?,
+            "ply" => ply_format::surface_mesh_from_ply(&input_file)?,
             _ => {
                 return Err(anyhow!(
                     "Unsupported file format extension \"{}\" for reading surface meshes",
diff --git a/splashsurf/src/io/ply_format.rs b/splashsurf/src/io/ply_format.rs
index 82679b8f..a2f43e70 100644
--- a/splashsurf/src/io/ply_format.rs
+++ b/splashsurf/src/io/ply_format.rs
@@ -4,6 +4,10 @@ use anyhow::{anyhow, Context};
 use ply_rs as ply;
 use ply_rs::ply::Property;
 
+use splashsurf_lib::mesh::AttributeData;
+use splashsurf_lib::mesh::MeshAttribute;
+use splashsurf_lib::mesh::MeshWithData;
+use splashsurf_lib::mesh::TriMesh3d;
 use splashsurf_lib::nalgebra::Vector3;
 use splashsurf_lib::Real;
 
@@ -49,3 +53,139 @@ pub fn particles_from_ply<R: Real, P: AsRef<Path>>(
 
     Ok(particles)
 }
+
+/// Tries to read a surface mesh from the VTK file at the given path
+pub fn surface_mesh_from_ply<R: Real, P: AsRef<Path>>(
+    ply_file: P,
+) -> Result<MeshWithData<R, TriMesh3d<R>>, anyhow::Error> {
+    let mut ply_file = std::fs::File::open(ply_file).unwrap();
+    let parser = ply::parser::Parser::<ply::ply::DefaultElement>::new();
+
+    let ply = parser
+        .read_ply(&mut ply_file)
+        .context("Failed to read PLY file")?;
+    let vertices_normals = ply
+        .payload
+        .get("vertex")
+        .ok_or(anyhow!("PLY file is missing a 'vertex' element"))?;
+
+    let vertices_normals: Vec<(Vector3<_>, Vector3<_>)> = vertices_normals
+        .into_iter()
+        .map(|e| {
+            let vertex = (
+                e.get("x").unwrap(),
+                e.get("y").unwrap(),
+                e.get("z").unwrap(),
+                e.get("nx").unwrap(),
+                e.get("ny").unwrap(),
+                e.get("nz").unwrap(),
+            );
+
+            let v = match vertex {
+                (
+                    Property::Float(x),
+                    Property::Float(y),
+                    Property::Float(z),
+                    Property::Float(nx),
+                    Property::Float(ny),
+                    Property::Float(nz),
+                ) => (
+                    Vector3::new(
+                        R::from_f32(*x).unwrap(),
+                        R::from_f32(*y).unwrap(),
+                        R::from_f32(*z).unwrap(),
+                    ),
+                    Vector3::new(
+                        R::from_f32(*nx).unwrap(),
+                        R::from_f32(*ny).unwrap(),
+                        R::from_f32(*nz).unwrap(),
+                    ),
+                ),
+                _ => {
+                    return Err(anyhow!(
+                        "Vertex properties have wrong PLY data type (expected float)"
+                    ))
+                }
+            };
+
+            Ok(v)
+        })
+        .map(|vn| vn.unwrap())
+        .collect();
+
+    let vertices: Vec<Vector3<_>> = vertices_normals.iter().map(|vn| vn.0.clone()).collect();
+    let normals: Vec<Vector3<_>> = vertices_normals.iter().map(|vn| vn.1.clone()).collect();
+
+    let faces = ply
+        .payload
+        .get("face")
+        .ok_or(anyhow!("PLY file is missing a 'face' element"))?;
+
+    let triangles = faces
+        .into_iter()
+        .map(|e| {
+            // This is as per what blender creates for a
+            let indices = e.get("vertex_indices");
+            if let Some(indices) = indices {
+                if let Property::ListUInt(indices) = indices {
+                    if indices.len() == 3 {
+                        return Ok([
+                            indices[0] as usize,
+                            indices[1] as usize,
+                            indices[2] as usize,
+                        ]);
+                    } else {
+                        return Err(anyhow!(
+                            "Invalid number of vertex indices per cell: {}",
+                            indices.len()
+                        ));
+                    }
+                } else {
+                    return Err(anyhow!(
+                        "Index properties have wrong PLY data type (expected uint)"
+                    ));
+                }
+            } else {
+                return Err(anyhow!(
+                    "Vertex properties have wrong PLY data type (expected uint)"
+                ));
+            }
+        })
+        .map(|e| e.unwrap())
+        .collect();
+
+    let normals = MeshAttribute::new("normals", AttributeData::Vector3Real(normals));
+    Ok(MeshWithData::new(TriMesh3d {
+        vertices,
+        triangles,
+    })
+    .with_point_data(normals))
+}
+
+#[cfg(test)]
+pub mod test {
+    use super::*;
+
+    #[test]
+    fn test_convert_cube() -> Result<(), anyhow::Error> {
+        let input_file = Path::new("../data/cube.ply");
+
+        let mesh: MeshWithData<f32, _> = surface_mesh_from_ply(input_file).with_context(|| {
+            format!(
+                "Failed to load surface mesh from file \"{}\"",
+                input_file.display()
+            )
+        })?;
+
+        assert_eq!(mesh.mesh.vertices.len(), 24);
+        assert_eq!(mesh.mesh.triangles.len(), 12);
+        let normals = mesh.point_attributes.iter().find(|a| a.name == "normals");
+        if let Some(MeshAttribute { data, .. }) = normals {
+            if let AttributeData::Vector3Real(normals) = data {
+                assert_eq!(normals.len(), 24)
+            }
+        }
+
+        Ok(())
+    }
+}