diff --git a/src/recipes.jl b/src/recipes.jl
index 952bf0aa..6901acb8 100644
--- a/src/recipes.jl
+++ b/src/recipes.jl
@@ -987,6 +987,74 @@ end
 end
 
 
+# ---------------------------------------------------------------------------
+# Histogram 3D
+
+# Check if all neighbours have the same value
+# TODO: Improve this somehow
+function _allneighbours(arr, I::CartesianIndex, val = 0.0)
+    arr[I] == val &&
+    arr[I + CartesianIndex(1,0)]  == val &&
+    arr[I + CartesianIndex(0,1)]  == val &&
+    arr[I - CartesianIndex(1,0)]  == val &&
+    arr[I - CartesianIndex(0,1)]  == val &&
+    arr[I + CartesianIndex(1,1)]  == val &&
+    arr[I + CartesianIndex(1,-1)] == val &&
+    arr[I - CartesianIndex(1,1)]  == val &&
+    arr[I - CartesianIndex(1,-1)] == val
+end
+
+@recipe function f(::Type{Val{:bins3d}}, x, y, z)
+    edge_x, edge_y, weights = x, y, z.surf
+
+    if !plotattributes[:show_empty_bins]
+        edge_x = repeat(edge_x,inner=(4))[2:end-1]
+        edge_y = repeat(edge_y,inner=(4))[2:end-1]
+        float_weights = spzeros(length(edge_x),length(edge_y))
+        float_weights[2:end-1,2:end-1] = repeat(float(weights),inner=(4,4))
+        m = falses(size(float_weights))
+        for I in CartesianIndices((2:length(edge_x)-1,2:length(edge_y)-1))
+            m[I] = _allneighbours(float_weights,I,0)
+        end
+        float_weights[m] .= NaN
+    else
+        edge_x = repeat(edge_x,inner=(2))
+        edge_y = repeat(edge_y,inner=(2))
+        float_weights = spzeros(length(edge_x),length(edge_y))
+        float_weights[2:end-1,2:end-1] .= repeat(float(weights),inner=(2,2))
+    end
+
+    x := edge_x
+    y := edge_y
+    z := Surface(permutedims(float_weights))
+    #seriescolor := cgrad([:blue,:blue])
+    colorbar := false
+    seriestype := :surface
+    linealpha := 1.0
+    linecolor := :black
+
+    ()
+
+end
+Plots.@deps bins3d surface wireframe
+
+
+@recipe function f(::Type{Val{:histogram3d}}, x, y, z)
+    h = _make_hist(
+        (x, y),
+        plotattributes[:bins],
+        normed = plotattributes[:normalize],
+        weights = plotattributes[:weights],
+    )
+    x := h.edges[1]
+    y := h.edges[2]
+    z := Surface(h.weights)
+    seriestype := :bins3d
+    ()
+end
+@deps histogram3d bins3d
+
+
 # ---------------------------------------------------------------------------
 # pie
 @recipe function f(::Type{Val{:pie}}, x, y, z)
@@ -1430,7 +1498,7 @@ end
 
 @recipe function f(shapes::AVec{Shape})
     seriestype --> :shape
-    # For backwards compatibility, column vectors of segmenting attributes are 
+    # For backwards compatibility, column vectors of segmenting attributes are
     # interpreted as having one element per shape
     for attr in union(_segmenting_array_attributes, _segmenting_vector_attributes)
         v = get(plotattributes, attr, nothing)
diff --git a/src/shorthands.jl b/src/shorthands.jl
index dfcc2302..b03e201a 100644
--- a/src/shorthands.jl
+++ b/src/shorthands.jl
@@ -107,6 +107,25 @@ julia> histogram2d(randn(10_000),randn(10_000))
 """
 @shorthands histogram2d
 
+"""
+    histogram3d(x,y)
+    histogram3d!(x,y)
+
+Plot a three-dimensional histogram.
+
+# Arguments
+
+- `bins`: Number of bins (if an `Integer`) or bin edges (if an `AbtractVector`)
+- `weights`: Vector of weights for the values in `x`. Each entry of x contributes
+             its weight to the height of its bin.
+
+# Example
+```julia-repl
+julia> histogram3d(randn(10_000),randn(10_000))
+```
+"""
+@shorthands histogram3d
+
 """
     density(x)
     density!(x)