diff --git a/src/Graphs.jl b/src/Graphs.jl
index c8d209cf5..18729063d 100644
--- a/src/Graphs.jl
+++ b/src/Graphs.jl
@@ -248,6 +248,7 @@ export
     transitivereduction,
     yen_k_shortest_paths,
     desopo_pape_shortest_paths,
+    bidijkstra_shortest_path,
 
     # centrality
     betweenness_centrality,
diff --git a/src/shortestpaths/dijkstra.jl b/src/shortestpaths/dijkstra.jl
index 39d58559b..e43d311c1 100644
--- a/src/shortestpaths/dijkstra.jl
+++ b/src/shortestpaths/dijkstra.jl
@@ -75,9 +75,8 @@ function dijkstra_shortest_paths(
     distmx::AbstractMatrix{T}=weights(g);
     allpaths=false,
     trackvertices=false,
-    maxdist=typemax(T)
-    ) where T <: Real where U <: Integer
-
+    maxdist=typemax(T),
+) where {T<:Real} where {U<:Integer}
     nvg = nv(g)
     dists = fill(typemax(T), nvg)
     parents = zeros(U, nvg)
@@ -163,9 +162,142 @@ function dijkstra_shortest_paths(
     distmx::AbstractMatrix=weights(g);
     allpaths=false,
     trackvertices=false,
-    maxdist=typemax(eltype(distmx))
+    maxdist=typemax(eltype(distmx)),
 )
     return dijkstra_shortest_paths(
         g, [src;], distmx; allpaths=allpaths, trackvertices=trackvertices, maxdist=maxdist
     )
 end
+
+function relax(u,
+               v, 
+               distmx::AbstractMatrix{T}, 
+               dists::Vector{T}, 
+               parents::Vector{U}, 
+               visited::Vector{Bool}, 
+               Q::PriorityQueue{U,T};
+               allpaths=false,
+               pathcounts=nothing,
+               preds=nothing,
+               forward=true
+) where {T<:Real} where {U<:Integer}
+    alt = dists[u] + (forward ? distmx[u, v] : distmx[v, u])
+
+    if !visited[v]
+        visited[v] = true
+        dists[v] = alt
+        parents[v] = u
+        
+        if !isnothing(pathcounts)
+            pathcounts[v] += pathcounts[u]
+        end
+        if allpaths
+            preds[v] = [u;]
+        end
+        Q[v] = alt
+    elseif alt < dists[v]
+        dists[v] = alt
+        parents[v] = u
+        #615
+        if !isnothing(pathcounts)
+            pathcounts[v] = pathcounts[u]
+        end
+        if allpaths
+            resize!(preds[v], 1)
+            preds[v][1] = u
+        end
+        Q[v] = alt
+    elseif alt == dists[v]
+        if !isnothing(pathcounts)
+            pathcounts[v] += pathcounts[u]
+        end
+        if allpaths
+            push!(preds[v], u)
+        end
+    end
+end
+
+"""
+    bidijkstra_shortest_paths(g, src, dst, distmx=weights(g));
+
+Perform [Bidirectional Dijkstra's algorithm](https://www.homepages.ucl.ac.uk/~ucahmto/math/2020/05/30/bidirectional-dijkstra.html)
+on a graph, computing the shortest path between `src` and `dst`.
+
+# Examples
+```jldoctest
+julia> using Graphs
+
+julia> bidijkstra_shortest_path(cycle_graph(5), 1, 4)
+3-element Vector{Int64}:
+ 1
+ 5
+ 4
+
+julia> bidijkstra_shortest_path(path_graph(5), 1, 4)
+4-element Vector{Int64}:
+ 1
+ 2
+ 3
+ 4
+```
+"""
+function bidijkstra_shortest_path(
+    g::AbstractGraph,
+    src::U,
+    dst::U,
+    distmx::AbstractMatrix{T}=weights(g)
+) where {T<:Real} where {U<:Integer}
+    if src == dst
+        return Int[]
+    end
+    # keep weight of the best seen path and the midpoint vertex
+    μ, mid_v = typemax(T), -1
+    nvg = nv(g)
+    dists_f, dists_b= fill(typemax(T), nvg), fill(typemax(T), nvg)
+    parents_f, parents_b= zeros(U, nvg), zeros(U, nvg)
+    visited_f, visited_b = zeros(Bool, nvg),zeros(Bool, nvg)
+    preds_f, preds_b = fill(Vector{U}(), nvg), fill(Vector{U}(), nvg)
+    Qf, Qb = PriorityQueue{U,T}(), PriorityQueue{U,T}()
+
+    dists_f[src], dists_b[dst]= zero(T), zero(T)
+    visited_f[src], visited_b[dst]= true, true
+    Qf[src], Qb[dst] = zero(T), zero(T)
+
+    while !isempty(Qf) && !isempty(Qb)
+        uf, ub = dequeue!(Qf), dequeue!(Qb)
+
+        for v in outneighbors(g, uf)
+            relax(uf, v, distmx, dists_f, parents_f, visited_f, Qf)
+            if visited_b[v] && (dists_f[uf]+distmx[uf,v]+dists_b[v]) < μ
+                # we have found an edge between the forward and backward exploration
+                μ = dists_f[uf]+distmx[uf,v]+dists_b[v]
+                mid_v = v
+            end
+        end
+
+        for v in inneighbors(g, ub)
+            relax(ub, v, distmx, dists_b, parents_b, visited_b, Qb; forward=false)
+            if visited_f[v] && (dists_f[v]+distmx[v,ub]+dists_b[ub]) < μ
+                # we have found an edge between the forward and backward exploration
+                μ = dists_f[v]+distmx[v,ub]+dists_b[ub]
+                mid_v = v
+            end
+        end
+        if dists_f[uf]+dists_b[ub] >= μ
+            break
+        end
+    end
+    if mid_v == -1
+        # no path exists between source and destination
+        return Int[]
+    end
+    ds_f = DijkstraState{T,U}(parents_f, dists_f, preds_f, zeros(nvg), Vector{U}())
+    ds_b = DijkstraState{T,U}(parents_b, dists_b, preds_b, zeros(nvg), Vector{U}())
+    if mid_v == src
+        return reverse(enumerate_paths(ds_b, mid_v))
+    elseif mid_v ==dst
+        return enumerate_paths(ds_f, mid_v)
+    end
+    return vcat(enumerate_paths(ds_f, mid_v), reverse(enumerate_paths(ds_b, mid_v)[1:end-1]))
+end
+
diff --git a/test/runtests.jl b/test/runtests.jl
index cbb8763bb..2ca2272c2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -105,6 +105,7 @@ tests = [
     "shortestpaths/bellman-ford",
     "shortestpaths/desopo-pape",
     "shortestpaths/dijkstra",
+    "shortestpaths/bidijkstra",
     "shortestpaths/johnson",
     "shortestpaths/floyd-warshall",
     "shortestpaths/yen",
diff --git a/test/shortestpaths/bidijkstra.jl b/test/shortestpaths/bidijkstra.jl
new file mode 100644
index 000000000..6c3005d6f
--- /dev/null
+++ b/test/shortestpaths/bidijkstra.jl
@@ -0,0 +1,19 @@
+@testset "Bidijkstra" begin
+    g3 = path_graph(5)
+    g4 = path_digraph(5)
+
+    d1 = float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0])
+    d2 = sparse(float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0]))
+    for g in testgraphs(g3), dg in testdigraphs(g4)
+        @test @inferred(bidijkstra_shortest_path(g, 1, 4, d1)) ==
+            @inferred(bidijkstra_shortest_path(dg, 1, 4, d1)) ==
+            @inferred(bidijkstra_shortest_path(g, 1, 4, d2))
+        @test isempty(@inferred(bidijkstra_shortest_path(dg, 4, 1)))
+    end
+
+    # test for #1258
+    g = complete_graph(4)
+    w = float([1 1 1 4; 1 1 1 1; 1 1 1 1; 4 1 1 1])
+    ds = dijkstra_shortest_paths(g, 1, w)
+    @test length(bidijkstra_shortest_path(g, 1, 4, w)) == 3 # path is a sequence of vertices
+end
diff --git a/test/shortestpaths/dijkstra.jl b/test/shortestpaths/dijkstra.jl
index a4f0f21b2..ec24c17a0 100644
--- a/test/shortestpaths/dijkstra.jl
+++ b/test/shortestpaths/dijkstra.jl
@@ -111,4 +111,20 @@
             ds = @inferred(dijkstra_shortest_paths(g, 3, m;maxdist=3.0))
             @test ds.dists == [2, 1, 0, Inf, Inf, 3]
         end
+
+    # bidijkstra_shortest_path
+    g4 = path_digraph(5)
+    d1 = float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0])
+    d2 = sparse(float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0]))
+
+    for g in testdigraphs(g4)
+        x = @inferred(dijkstra_shortest_paths(g, 2, d1))
+        p = enumerate_paths(x, 4)
+        y = @inferred(bidijkstra_shortest_path(g, 2, 4, d1))
+        z = @inferred(bidijkstra_shortest_path(g, 2, 4, d2))
+        
+        @test p == y == z
+    end
+
+   
 end