diff --git a/src/args.jl b/src/args.jl
index 5e624c0c..7ea061e2 100644
--- a/src/args.jl
+++ b/src/args.jl
@@ -51,12 +51,12 @@ const _axesAliases = Dict{Symbol,Symbol}(
 )
 
 const _3dTypes = [
-    :path3d, :scatter3d, :surface, :wireframe, :contour3d, :volume, :mesh3d
+    :path3d, :scatter3d, :surface, :wireframe, :contour3d, :volume, :mesh3d, :histogram3d
 ]
 const _allTypes = vcat([
     :none, :line, :path, :steppre, :stepmid, :steppost, :sticks, :scatter,
     :heatmap, :hexbin, :barbins, :barhist, :histogram, :scatterbins,
-    :scatterhist, :stepbins, :stephist, :bins2d, :histogram2d, :histogram3d,
+    :scatterhist, :stepbins, :stephist, :bins2d, :histogram2d,
     :density, :bar, :hline, :vline,
     :contour, :pie, :shape, :image
 ], _3dTypes)
@@ -94,13 +94,15 @@ const _typeAliases = Dict{Symbol,Symbol}(
     :imagesc       => :image,
     :hist          => :histogram,
     :hist2d        => :histogram2d,
+    :hist3d        => :histogram3d,
+    :lego          => :histogram3d,
     :bezier        => :curves,
     :bezier_curves => :curves,
 )
 
 add_non_underscore_aliases!(_typeAliases)
 
-const _histogram_like = [:histogram, :barhist, :barbins]
+const _histogram_like = [:histogram, :barhist, :barbins, :histogram3d]
 const _line_like = [:line, :path, :steppre, :stepmid, :steppost]
 const _surface_like = [:contour, :contourf, :contour3d, :heatmap, :surface, :wireframe, :image]
 
diff --git a/src/axes.jl b/src/axes.jl
index fc99e5b8..caf51d9b 100644
--- a/src/axes.jl
+++ b/src/axes.jl
@@ -521,7 +521,7 @@ function widen(lmin, lmax, scale=:identity)
 end
 
 # figure out if widening is a good idea.
-const _widen_seriestypes = (:line, :path, :steppre, :stepmid, :steppost, :sticks, :scatter, :barbins, :barhist, :histogram, :scatterbins, :scatterhist, :stepbins, :stephist, :bins2d, :histogram2d, :bar, :shape, :path3d, :scatter3d)
+const _widen_seriestypes = (:line, :path, :steppre, :stepmid, :steppost, :sticks, :scatter, :barbins, :barhist, :histogram, :scatterbins, :scatterhist, :stepbins, :stephist, :bins2d, :histogram2d, :bar, :shape, :path3d, :scatter3d, :histogram3d)
 
 function default_should_widen(axis::Axis)
     if axis[:widen] isa Bool
diff --git a/src/recipes.jl b/src/recipes.jl
index c2911db9..c2213dd2 100644
--- a/src/recipes.jl
+++ b/src/recipes.jl
@@ -987,6 +987,118 @@ 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
+
+    # Create the bins
+    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
+
+    # pyplot can't handle sparse arrays it seems
+    # Maybe get rid off sparse arrays at all, but if any bins-entry is large the additional points for plotting may use up a lot of memory otherwise for zeros and NaNs
+    if backend_name() == :pyplot
+        float_weights = Surface(permutedims(Array(float_weights)))
+    else
+        float_weights = Surface(permutedims(float_weights))
+    end
+
+    # Handle input of one seriescolor only (maybe exchangeable by get_series_color()?)
+    if isa(plotattributes[:seriescolor], Symbol) && plotattributes[:seriescolor] != :auto
+        seriescolor := cgrad([plotattributes[:seriescolor], plotattributes[:seriescolor]])
+    end
+
+    seriestype := :surface
+    colorbar := false
+    x := edge_x
+    y := edge_y
+    z := float_weights
+
+    ()
+end
+Plots.@deps bins3d surface #wireframe mesh3d
+
+
+
+@recipe function f(::Type{Val{:bricks3d}}, x, y, z)
+    edge_x, edge_y, weights = x, y, z.surf
+
+    x_step = edge_x[2] - edge_x[1]
+    y_step = edge_y[2] - edge_y[1]
+    x_len = length(x)
+    y_len = length(y)
+
+    temp_x = vec([0.0 0.0 1.0 1.0]).*x_step .+ edge_x'
+    temp_y = vec([1.0 1.0 0.0 0.0]).*y_step .+ edge_y'
+
+    z_help =     [1.0 0.0 0.0 1.0;
+                  0.0 1.0 1.0 0.0;
+                  0.0 1.0 1.0 0.0;
+                  1.0 0.0 0.0 1.0];
+
+    for (i, c) in enumerate(weights)
+        itx = (i-1)%y_len + 1
+        ity = floor(Int64,i/x_len) + 1
+        @series begin
+            seriestype := :surface
+            x := temp_x[:,itx]
+            y := temp_y[:,ity]
+            z := c == 0 ? Surface(permutedims(NaN .* z_help)) : Surface(permutedims(c .* z_help))
+            ()
+        end
+    end
+
+
+end
+Plots.@deps bins3d surface #wireframe mesh3d
+
+
+@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 := :bricks3d#:bins3d
+    ()
+end
+@deps histogram3d bins3d bricks3d
+
+
 # ---------------------------------------------------------------------------
 # pie
 @recipe function f(::Type{Val{:pie}}, x, y, z)
diff --git a/src/shorthands.jl b/src/shorthands.jl
index 74a7b26c..36c84590 100644
--- a/src/shorthands.jl
+++ b/src/shorthands.jl
@@ -107,6 +107,26 @@ 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
+@shorthands legoplot
+
 """
     density(x)
     density!(x)
diff --git a/src/subplots.jl b/src/subplots.jl
index 35167349..ae893a85 100644
--- a/src/subplots.jl
+++ b/src/subplots.jl
@@ -49,6 +49,7 @@ function should_add_to_legend(series::Series)
             :hexbin,
             :bins2d,
             :histogram2d,
+            :histogram3d,
             :hline,
             :vline,
             :contour,