diff --git a/src/MeshIO.jl b/src/MeshIO.jl
index 01d0c24..0f50e29 100644
--- a/src/MeshIO.jl
+++ b/src/MeshIO.jl
@@ -9,6 +9,8 @@ using FileIO: FileIO, @format_str, Stream, File, stream, skipmagic
 
 import Base.show
 
+include("util.jl")
+
 include("io/off.jl")
 include("io/ply.jl")
 include("io/stl.jl")
diff --git a/src/io/obj.jl b/src/io/obj.jl
index 0ebf193..5f48873 100644
--- a/src/io/obj.jl
+++ b/src/io/obj.jl
@@ -61,70 +61,34 @@ function load(io::Stream{format"OBJ"}; facetype=GLTriangleFace,
     point_attributes = Dict{Symbol, Any}()
     non_empty_faces = filtertuple(!isempty, f_uv_n_faces)
 
+    # Do we have faces with different indices for positions and normals 
+    # (and texture coordinates) per vertex?
     if length(non_empty_faces) > 1
-        N = length(points)
-        void = tuple((_typemax(eltype(facetype)) for _ in 1:length(non_empty_faces))...)
-        vertices = fill(void, N)
 
-        if !isempty(v_normals)
-            point_attributes[:normals] = Vector{normaltype}(undef, N)
-        end
+        # map vertices with distinct indices for possition and normal (and uv)
+        # to new indices, updating faces along the way
+        faces, attrib_maps = merge_vertex_attribute_indices(non_empty_faces)
+
+        # Update order of vertex attributes
+        points = points[attrib_maps[1]]
+        counter = 2
         if !isempty(uv)
-            point_attributes[:uv] = Vector{uvtype}(undef, N)
+            point_attributes[:uv] = uv[attrib_maps[counter]]
+            counter += 1
         end
-
-        for (k, fs) in enumerate(zip(non_empty_faces...))
-            f = collect(first(fs)) # position indices
-            for i in eachindex(non_empty_faces)
-                l = 2
-                vertex = getindex.(fs, i) # one of each indices (pos/uv/normal)
-
-                if vertices[vertex[1]] == void
-                    # Replace void
-                    vertices[vertex[1]] = vertex
-                    f[i] = vertex[1]
-                    if !isempty(uv)
-                        point_attributes[:uv][vertex[1]] = uv[vertex[l]]
-                        l += 1
-                    end
-                    if !isempty(v_normals)
-                        point_attributes[:normals][vertex[1]] = v_normals[vertex[l]]
-                    end
-                elseif vertices[vertex[1]] == vertex
-                    # vertex is correct, nothing to replace
-                    f[i] = vertex[1]
-                else
-                    @views j = findfirst(==(vertex), vertices[N+1:end])
-                    if j === nothing
-                        # vertex is unique, add it as a new one and adjust
-                        # points, uv, normals
-                        push!(vertices, vertex)
-                        f[i] = length(vertices)
-                        push!(points, points[vertex[1]])
-                        if !isempty(uv)
-                            push!(point_attributes[:uv], uv[vertex[l]])
-                            l += 1
-                        end
-                        if !isempty(v_normals)
-                            push!(point_attributes[:normals], v_normals[vertex[l]])
-                        end
-                    else
-                        # vertex has already been added, adjust face
-                        # (points, uv, normals correct because they've been pushed)
-                        f[i] = j + N
-                    end
-                end
-            end
-            # remap indices
-            faces[k] = facetype(f)
+        if !isempty(v_normals)
+            point_attributes[:normals] = v_normals[attrib_maps[counter]]
         end
+
     else # we have vertex indexing - no need to remap
+
         if !isempty(v_normals)
             point_attributes[:normals] = v_normals
         end
         if !isempty(uv)
             point_attributes[:uv] = uv
         end
+        
     end
 
     return Mesh(meta(points; point_attributes...), faces)
diff --git a/src/util.jl b/src/util.jl
new file mode 100644
index 0000000..1c24a8f
--- /dev/null
+++ b/src/util.jl
@@ -0,0 +1,56 @@
+# Graphics backends like OpenGL only have one index buffer so the indices to
+# positions, normals and texture coordinates cannot be different. E.g. a face
+# cannot use positional indices (1, 2, 3) and normal indices (1, 1, 2). In that
+# case we need to remap normals such that new_normals[1, 2, 3] = normals[[1, 1, 2]]
+
+
+# ...
+_typemin(x) = typemin(x)
+_typemin(::Type{OffsetInteger{N, T}}) where {N, T} = typemin(T) - N
+
+merge_vertex_attribute_indices(faces...) = merge_vertex_attribute_indices(faces)
+
+function merge_vertex_attribute_indices(faces::Tuple)
+    FaceType = eltype(faces[1])
+    IndexType = eltype(FaceType)
+    D = length(faces)
+    N = length(faces[1])
+
+    # (pos_idx, normal_idx, uv_idx, ...) -> new_idx
+    vertex_index_map = Dict{NTuple{D, UInt32}, IndexType}()
+    # faces after remapping (0 based assumed)
+    new_faces = sizehint!(FaceType[], N)
+    temp = IndexType[]           # keeping track of vertex indices of a face
+    counter = _typemin(IndexType)
+    # for remaping attribs, i.e. `new_attrib = old_attrib[index2vertex[attrib_index]]`
+    index2vertex = ntuple(_ -> sizehint!(UInt32[], N), D)
+
+    for i in eachindex(faces[1])
+        # (pos_faces[i], normal_faces[i], uv_faces[i], ...)
+        attrib_faces = getindex.(faces, i)
+        empty!(temp)
+        
+        for j in eachindex(attrib_faces[1])
+            # (pos_index, normal_idx, uv_idx, ...)
+            # = (pos_faces[i][j], normal_faces[i][j], uv_faces[i][j], ...)
+            vertex = GeometryBasics.value.(getindex.(attrib_faces, j)) # 1 based
+
+            # if combination of indices in vertex is new, make a new index
+            if !haskey(vertex_index_map, vertex)
+                vertex_index_map[vertex] = counter
+                counter = IndexType(counter + 1)
+                push!.(index2vertex, vertex)
+            end
+
+            # keep track of the (new) index for this vertex
+            push!(temp, vertex_index_map[vertex])
+        end
+
+        # store face with new indices
+        push!(new_faces, FaceType(temp...))
+    end
+
+    sizehint!(new_faces, length(new_faces))
+
+    return new_faces, index2vertex
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index d9a7edd..70ae35b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -133,8 +133,8 @@ end
         @testset "OBJ" begin
             msh = load(joinpath(tf, "test.obj"))
             @test length(faces(msh)) == 3954
-            @test length(coordinates(msh)) == 2520
-            @test length(normals(msh)) == 2520
+            @test length(coordinates(msh)) == 2519
+            @test length(normals(msh)) == 2519
             @test test_face_indices(msh)
 
             msh = load(joinpath(tf, "cube.obj")) # quads
@@ -176,5 +176,25 @@ end
             #@test typeof(msh) == GLNormalMesh
             #test_face_indices(msh)
         end
+
+        @testset "Index remapping" begin
+            pos_faces    = GLTriangleFace[(5, 6, 7), (5, 6, 8), (5, 7, 8)]
+            normal_faces = GLTriangleFace[(5, 6, 7), (3, 6, 8), (5, 7, 8)]
+            uv_faces     = GLTriangleFace[(1, 2, 3), (4, 2, 5), (1, 3, 1)]
+            
+            #   unique combinations    ->      new indices
+            # 551 662 773 534 885 881      1 2 3 4 5 6 (or 0..5 with 0 based indices)
+            faces, maps = MeshIO.merge_vertex_attribute_indices(pos_faces, normal_faces, uv_faces)
+
+            @test length(faces) == 3
+            @test faces == GLTriangleFace[(1, 2, 3), (4, 2, 5), (1, 3, 6)]
+
+            # maps are structured as map[new_index] = old_index, so they grab the
+            # first/second/third index of the unique combinations above
+            # maps = (pos_map, normal_map, uv_map)
+            @test maps[1] == [5, 6, 7, 5, 8, 8]
+            @test maps[2] == [5, 6, 7, 3, 8, 8]
+            @test maps[3] == [1, 2, 3, 4, 5, 1]
+        end
     end
 end