diff --git a/src/models.c b/src/models.c
index e34099e6..55621ee1 100644
--- a/src/models.c
+++ b/src/models.c
@@ -2473,6 +2473,74 @@ void MeshBinormals(Mesh *mesh)
     }
 }
 
+// Smooth (average) vertex normals
+void MeshNormalsSmooth(Mesh *mesh)
+{
+    #define EPSILON 0.000001 // A small number
+
+    int uvCounter = 0;
+    Vector3 *uniqueVertices = (Vector3 *)RL_CALLOC(mesh->vertexCount, sizeof(Vector3));
+    Vector3 *summedNormals = (Vector3 *)RL_CALLOC(mesh->vertexCount, sizeof(Vector3));
+
+    int uiCounter = 0;
+    int *uniqueIndices = (int *)RL_CALLOC(mesh->vertexCount, sizeof(int));
+
+    // Sum normals grouped by vertex
+    for (int i = 0; i < mesh->vertexCount; i++)
+    {
+        Vector3 v = { mesh->vertices[(i + 0)*3 + 0], mesh->vertices[(i + 0)*3 + 1], mesh->vertices[(i + 0)*3 + 2] };
+        Vector3 n = { mesh->normals[(i + 0)*3 + 0], mesh->normals[(i + 0)*3 + 1], mesh->normals[(i + 0)*3 + 2] };
+
+        bool matched = false;
+
+        // TODO: Matching vertices is brute force O(N^2). Do it more efficiently?
+        for (int j = 0; j < uvCounter; j++)
+        {
+            Vector3 uv = uniqueVertices[j];
+
+            bool match = true;
+            match = match && fabs(uv.x - v.x) < EPSILON;
+            match = match && fabs(uv.y - v.y) < EPSILON;
+            match = match && fabs(uv.z - v.z) < EPSILON;
+
+            if (match)
+            {
+                matched = true;
+                summedNormals[j] = Vector3Add(summedNormals[j], n);
+                uniqueIndices[i] = j;
+                break;
+            }
+        }
+
+        if (!matched)
+        {
+            int j = uvCounter++;
+            uniqueVertices[j] = v;
+            summedNormals[j] = n;
+            uniqueIndices[i] = j;
+        }
+    }
+
+    // Average and update normals
+    for (int i = 0; i < mesh->vertexCount; i++)
+    {
+        int j = uniqueIndices[i];
+        Vector3 n = Vector3Normalize(summedNormals[j]);
+        mesh->normals[(i + 0)*3 + 0] = n.x;
+        mesh->normals[(i + 0)*3 + 1] = n.y;
+        mesh->normals[(i + 0)*3 + 2] = n.z;
+    }
+
+    // 2=normals, see rlUpdateMeshAt()
+    rlUpdateMesh(*mesh, 2, mesh->vertexCount);
+
+    RL_FREE(uniqueVertices);
+    RL_FREE(summedNormals);
+    RL_FREE(uniqueIndices);
+
+    TRACELOG(LOG_INFO, "MESH: Normals smoothed (%d vertices, %d unique)", mesh->vertexCount, uvCounter);
+}
+
 // Draw a model (with texture if set)
 void DrawModel(Model model, Vector3 position, float scale, Color tint)
 {
diff --git a/src/raylib.h b/src/raylib.h
index d38d53a0..2329305e 100644
--- a/src/raylib.h
+++ b/src/raylib.h
@@ -1325,6 +1325,7 @@ RLAPI Mesh GenMeshCubicmap(Image cubicmap, Vector3 cubeSize);
 RLAPI BoundingBox MeshBoundingBox(Mesh mesh);                                                           // Compute mesh bounding box limits
 RLAPI void MeshTangents(Mesh *mesh);                                                                    // Compute mesh tangents
 RLAPI void MeshBinormals(Mesh *mesh);                                                                   // Compute mesh binormals
+RLAPI void MeshNormalsSmooth(Mesh *mesh);                                                               // Smooth (average) vertex normals
 
 // Model drawing functions
 RLAPI void DrawModel(Model model, Vector3 position, float scale, Color tint);                           // Draw a model (with texture if set)