diff --git a/ext/PlottingExt.jl b/ext/PlottingExt.jl
index f63f43c57..111c0e5b8 100644
--- a/ext/PlottingExt.jl
+++ b/ext/PlottingExt.jl
@@ -1175,6 +1175,58 @@ function Sunny.plot_intensities!(panel, res::Sunny.PowderIntensities{Float64}; c
     return ax
 end
 
+function Sunny.plot_intensities!(panel, res::Sunny.BinnedIntensities{Float64}; colormap=nothing, colorrange=nothing, saturation=0.9, allpositive=true, units=nothing, into=nothing, title="", axisopts=NamedTuple(),divide_counts = true)
+
+    axisopts = (; title, axisopts...)
+    unit_energy, elabel = get_unit_energy(units, into)
+
+    data = if divide_counts
+      res.data ./ res.counts
+    else
+      res.data
+    end
+
+    data_nonan = copy(data)
+    good_dats = data_nonan[isfinite.(data_nonan)]
+    data_nonan[.!isfinite.(data_nonan)] .= sum(good_dats) / length(good_dats)
+
+    colorrange_suggest = colorrange_from_data(; data=reshape(data_nonan, 1, size(data)...), saturation, sensitivity=0, allpositive)
+    colormap = @something colormap (allpositive ? :gnuplot2 : :bwr)
+    colorrange = @something colorrange colorrange_suggest
+
+    # Low-dimension cases
+    n_dims_resolved = count(res.params.numbins .!= 1)
+
+    ax = if n_dims_resolved == 0
+      # No resolved data: Just display the one value!
+      ax = Makie.Axis(panel[1,1];axisopts...)
+      Makie.text!(ax,0,0;text = "$(Sunny.number_to_simple_string(data[1];digits = 4))")
+      ax
+    elseif n_dims_resolved == 1
+      # Only resolved on one axis!
+      x_axis = findfirst(res.params.numbins .!= 1)
+      ax = Makie.Axis(panel[1,1];xlabel = x_axis == 4 ? elabel : Sunny.covector_name(res.params.covectors[x_axis,:]),ylabel = "Integrated Intensity",axisopts...)
+      bcs = Sunny.axes_bincenters(res.params)
+      bcs[4] /= unit_energy
+      Makie.barplot!(ax,bcs[x_axis],data[:];colormap, colorrange)
+      ax
+    elseif n_dims_resolved == 2
+      x_axis = findfirst(res.params.numbins .!= 1)
+      y_axis = x_axis + findfirst(res.params.numbins[x_axis+1:end] .!= 1)
+      ax = Makie.Axis(panel[1,1];
+                      xlabel = Sunny.covector_name(res.params.covectors[x_axis,:]),
+                      ylabel = y_axis == 4 ? elabel : Sunny.covector_name(res.params.covectors[y_axis,:]),axisopts...)
+      bcs = Sunny.axes_bincenters(res.params)
+      bcs[4] /= unit_energy
+      hm = Makie.heatmap!(ax,bcs[x_axis],bcs[y_axis],reshape(data,size(data,x_axis),size(data,y_axis));colormap, colorrange)
+      Makie.Colorbar(panel[1, 2], hm)
+      ax
+    elseif n_dims_resolved > 2
+      error("Higher-dimensional binned data visualization not yet implemented!")
+    end
+    ax
+end
+
 #=
   * `axisopts`: An additional collection of named arguments that will be passed
     to the `Makie.Axis` constructor. This allows to override the axis `xlabel`,
diff --git a/src/Binning/Binning.jl b/src/Binning/Binning.jl
index 6f76b377c..53ccfeff3 100644
--- a/src/Binning/Binning.jl
+++ b/src/Binning/Binning.jl
@@ -49,148 +49,150 @@ function Base.show(io::IO, ::MIME"text/plain", params::BinningParameters)
         first_edges = map(x -> x[1],bin_edges)
         last_edges = map(x -> x[end],bin_edges)
         @printf(io," from %+.3f to %+.3f along [", first_edges[k], last_edges[k])
-        axes_names = ["x","y","z","E"]
-        inMiddle = false
-        for j = 1:4
-            if params.covectors[k,j] != 0.
-                if(inMiddle)
-                    print(io," ")
-                end
-                @printf(io,"%+.2f d%s",params.covectors[k,j],axes_names[j])
-                inMiddle = true
-            end
-        end
+        print(io,covector_name(params.covectors[k,:]))
         @printf(io,"] (Δ = %.3f)", params.binwidth[k]/norm(params.covectors[k,:]))
         println(io)
     end
 end
 
-#=
-
-Base.copy(p::BinningParameters) = BinningParameters(copy(p.binstart),copy(p.binend),copy(p.binwidth),copy(p.covectors))
-
-# Support numbins as a (virtual) property, even though only the binwidth is stored
-Base.getproperty(params::BinningParameters, sym::Symbol) = sym == :numbins ? count_bins(params.binstart,params.binend,params.binwidth) : getfield(params,sym)
-
-function Base.setproperty!(params::BinningParameters, sym::Symbol, numbins)
-    if sym == :numbins
-        # *Ensure* that the last bin contains params.binend
-        params.binwidth .= (params.binend .- params.binstart) ./ (numbins .- 0.5)
-    else
-        setfield!(params,sym,numbins)
-    end
+function covector_name(cov)
+  str = ""
+  axes_names = ["x","y","z","E"]
+  inMiddle = false
+  for j = 1:4
+      if cov[j] != 0.
+          if(inMiddle)
+              str *= " "
+          end
+          str *= @sprintf("%+.2f d%s",cov[j],axes_names[j])
+          inMiddle = true
+      end
+  end
+  str
 end
 
-"""
-    count_bins(binstart,binend,binwidth)
+# Creates a binning scheme centered on the q_path, with the specified transverse
+# binning directions and bin widths.
+function specify_transverse_binning(q_path::QPath, first_transverse_axis, second_transverse_axis, first_transverse_binwidth, second_transverse_binwidth)
+  # Ensure path is non-empty and single-segment
+  if length(q_path.qs) < 2
+    error("Q path must have at least two points to infer bin width")
+  end
 
-Returns the number of bins in the binning scheme implied by `binstart`, `binend`, and `binwidth`.
-To count the bins in a [`BinningParameters`](@ref), use `params.numbins`.
+  Δq = q_path.qs[2] - q_path.qs[1]
 
-This function defines how partial bins are handled, so it should be used preferentially over
-computing the number of bins manually.
-"""
-count_bins(bin_start,bin_end,bin_width) = ceil.(Int64,(bin_end .- bin_start) ./ bin_width)
+  if !all([Δq ≈ dq for dq = diff(q_path.qs)])
+    error("Q path is multi-segment or irregular!")
+  end
 
-function BinningParameters(binstart,binend,binwidth;covectors = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1])
-    return BinningParameters(binstart,binend,binwidth,covectors)
-end
+  binstart = zero(MVector{4,Float64})
+  binstart[4] = -Inf # Default to integrate over all energies
 
-function BinningParameters(binstart,binend;numbins,kwargs...)
-    params = BinningParameters(binstart,binend,[0.,0,0,0];kwargs...)
-    params.numbins = numbins # Use the setproperty to do it correctly
-    params
-end
+  binend = zero(MVector{4,Float64})
+  binend[4] = 0
 
-"""
-    integrate_axes!(params::BinningParameters; axes)
-Integrate over one or more axes of the histogram by setting the number of bins
-in that axis to 1. Examples:
+  covectors = zero(MMatrix{4,4,Float64})
+  recip_directions = zeros(Float64,3,3)
+  recip_directions[:,1] .= Δq ./ norm(Δq)
+  recip_directions[:,2] .= first_transverse_axis
+  recip_directions[:,3] .= second_transverse_axis
 
-    integrate_axes!(params; axes = [2,3])
-    integrate_axes!(params; axes = 2)
-"""
-function integrate_axes!(params::BinningParameters;axes)
-    for k in axes
-        nbins = [params.numbins.data...]
-        nbins[k] = 1
-        params.numbins = SVector{4}(nbins)
-    end
-    return params
-end
+  if minimum(svd(recip_directions).S) < 1e-8
+    error("Axes are collinear!")
+  end
+  covectors[1:3,1:3] .= inv(recip_directions)
 
-# Find an axis-aligned bounding box containing the histogram
-function binning_parameters_aabb(params)
-    (; binstart, binend, covectors) = params
-    bin_edges = axes_binedges(params)
-    first_edges = map(x -> x[1],bin_edges)
-    last_edges = map(x -> x[end],bin_edges)
-    bin_edges = [first_edges last_edges]
-    this_corner = MVector{4,Float64}(undef)
-    q_corners = MMatrix{4,16,Float64}(undef)
-    for j = 1:16 # The sixteen corners of a 4-cube
-        for k = 1:4 # The four axes
-            this_corner[k] = bin_edges[k,1 + (j >> (k-1) & 1)]
-        end
-        q_corners[:,j] = covectors \ this_corner
-    end
-    lower_aabb_q = minimum(q_corners,dims=2)[1:3]
-    upper_aabb_q = maximum(q_corners,dims=2)[1:3]
-    return lower_aabb_q, upper_aabb_q
-end
-
-"""
-If `params` expects to bin values `(k,ω)` in absolute units, then calling
+  coords_start = covectors[1:3,1:3] * q_path.qs[1]
+  coords_end = covectors[1:3,1:3] * q_path.qs[end]
 
-    bin_rlu_as_absolute_units!(params::BinningParameters,[reciprocal lattice vectors])
+  binwidth = zero(MVector{4,Float64})
+  binwidth[1] = (covectors[1:3,1:3] * Δq)[1]
+  binwidth[2] = first_transverse_binwidth
+  binwidth[3] = second_transverse_binwidth
+  binwidth[4] = Inf
 
-will modifiy the `covectors` in `params` so that they will accept `(q,ω)` in Reciprocal Lattice Units (R.L.U.) instead.
-Conversly, if `params` expects `(q,ω)` R.L.U., calling
+  binstart[1:3] .= coords_start[1:3] .- binwidth[1:3]/2
+  binend[1:3] .= coords_end[1:3]
 
-    bin_absolute_units_as_rlu!(params::BinningParameters,[reciprocal lattice vectors])
+  # Check the original Q points end up in bin centers
+  in_bin(q) = (covectors[1:3,1:3] * q .- binstart[1:3]) ./ binwidth[1:3]
+  centering_error(q) = (in_bin(q) .- 0.5) .- (round.(Int64,in_bin(q) .- 0.5))
+  @assert all([norm(centering_error(q_path.qs[i])) < 1e-12 for i = 1:length(q_path.qs)])
 
-will adjust `params` to instead accept `(k,ω)` absolute units.
-
-The second argument may be a 3x3 matrix specifying the reciprocal lattice vectors, or a [`Crystal`](@ref).
-"""
-bin_absolute_units_as_rlu!, bin_rlu_as_absolute_units!
+  # Energy axis
+  covectors[4,:] .= [0,0,0,1]
 
-function bin_rlu_as_absolute_units!(params::BinningParameters,recip_vecs::AbstractMatrix)
-    covectorsK = params.covectors
-
-    # covectorsQ * q = covectorsK *  recip_vecs * q = covectorsK * k
-    # covectorsQ     = covectorsK *  recip_vecs
-    covectorsQ       = covectorsK * [recip_vecs [0;0;0]; [0 0 0] 1]
-    params.covectors = MMatrix{4,4}(covectorsQ)
-    params
+  BinningParameters(binstart,binend,binwidth,covectors)
 end
 
-# covectorsK * k = covectorsQ * inv(recip_vecs) * k = covectorsQ * q
-# covectorsK     = covectorsQ * inv(recip_vecs)
-bin_absolute_units_as_rlu!(params::BinningParameters,recip_vecs::AbstractMatrix) = bin_rlu_as_absolute_units!(params,inv(recip_vecs))
-
-bin_absolute_units_as_rlu!(params::BinningParameters,crystal::Crystal) = bin_absolute_units_as_rlu!(params,2π*inv(crystal.latvecs)')
-
-bin_rlu_as_absolute_units!(params::BinningParameters,crystal::Crystal) = bin_absolute_units_as_rlu!(params,crystal.latvecs'/2π)
+# Creates a binning scheme centered on the q_grid, with the specified transverse
+# binning direction and bin width.
+function specify_transverse_binning(q_grid::QGrid{2}, transverse_axis, transverse_binwidth)
+  # Ensure grid is non-empty and single-segment
+  if size(q_grid.qs,1) < 2 || size(q_grid.qs,2) < 2
+    error("2D Q grid must have at least two points in each direction")
+  end
+
+  Δq1 = q_grid.qs[2,1] - q_grid.qs[1,1]
+  Δq2 = q_grid.qs[1,2] - q_grid.qs[1,1]
+
+  if !all([Δq1 ≈ dq for dq = diff(q_grid.qs,dims=1)]) || !all([Δq2 ≈ dq for dq = diff(q_grid.qs,dims = 2)])
+    error("2D Q grid is irregular!")
+  end
+
+  binstart = zero(MVector{4,Float64})
+  binstart[4] = -Inf # Default to integrate over all energies
+
+  binend = zero(MVector{4,Float64})
+  binend[4] = 0
+
+  covectors = zero(MMatrix{4,4,Float64})
+  recip_directions = zeros(Float64,3,3)
+  recip_directions[:,1] .= Δq1 ./ norm(Δq1)
+  recip_directions[:,2] .= Δq2 ./ norm(Δq2)
+  recip_directions[:,3] .= transverse_axis
+
+  if minimum(svd(recip_directions).S) < 1e-8
+    error("Transverse axis is in-plane!")
+  end
+  covectors[1:3,1:3] .= inv(recip_directions)
+
+  coords_start = covectors[1:3,1:3] * q_grid.qs[1]
+  coords_end = covectors[1:3,1:3] * q_grid.qs[end]
+
+  first_binwidth = (covectors[1:3,1:3] * Δq1)[1]
+  second_binwidth = (covectors[1:3,1:3] * Δq2)[2]
+
+  binwidth = zero(MVector{4,Float64})
+  binwidth[1] = first_binwidth
+  binwidth[2] = second_binwidth
+  binwidth[3] = transverse_binwidth
+  binwidth[4] = Inf
+
+  binstart[1:3] .= coords_start[1:3] .- binwidth[1:3]/2
+  binend[1:3] .= coords_end[1:3]
+
+  # Check the original Q points end up in bin centers
+  in_bin(q) = (covectors[1:3,1:3] * q .- binstart[1:3]) ./ binwidth[1:3]
+  centering_error(q) = (in_bin(q) .- 0.5) .- (round.(Int64,in_bin(q) .- 0.5))
+  @assert all([norm(centering_error(q_grid.qs[i])) < 1e-8 for i = 1:length(q_grid.qs)])
+
+  # Energy axis
+  covectors[4,:] .= [0,0,0,1]
+  BinningParameters(binstart,binend,binwidth,covectors)
+end
 
 """
     unit_resolution_binning_parameters(sc::SampledCorrelations)
 
 Create [`BinningParameters`](@ref) which place one histogram bin centered at each possible `(q,ω)` scattering vector of the crystal.
 This is the finest possible binning without creating bins with zero scattering vectors in them.
-
-This function can be used without reference to a [`SampledCorrelations`](@ref) using an alternate syntax to manually specify the bin centers for the energy axis and the lattice size:
-
-    unit_resolution_binning_parameters(ω_bincenters,dims,[reciprocal lattice vectors])
-
-The last argument may be a 3x3 matrix specifying the reciprocal lattice vectors, or a [`Crystal`](@ref).
-
-Lastly, binning parameters for a single axis may be specifed by their bin centers:
-
-    (binstart,binend,binwidth) = unit_resolution_binning_parameters(bincenters::Vector{Float64})
 """
-function unit_resolution_binning_parameters(ωvals,dims,args...)
-    numbins = (dims...,length(ωvals))
+function unit_resolution_binning_parameters(sc;negative_energies = true)
+    ωvals = available_energies_including_zero(sc;negative_energies)
+
+    good_qs = available_wave_vectors(sc)
+    numbins = (size(good_qs)...,length(ωvals))
     # Bin centers should be at Sunny scattering vectors
     maxQ = 1 .- (1 ./ numbins)
     
@@ -202,7 +204,7 @@ function unit_resolution_binning_parameters(ωvals,dims,args...)
     binstart = (0.,0.,0.,minimum(ωvals)) .- (binwidth ./ 2)
     binend = (maxQ[1],maxQ[2],maxQ[3],maximum(ωvals)) # bin end is well inside of last bin
 
-    params = BinningParameters(binstart,binend,binwidth)
+    params = BinningParameters(binstart,binend,binwidth,I(4))
 
     # Special case for when there is only one bin in a direction
     for i = 1:4
@@ -215,81 +217,35 @@ function unit_resolution_binning_parameters(ωvals,dims,args...)
     params
 end
 
-unit_resolution_binning_parameters(sc::SampledCorrelations; negative_energies=false,kwargs...) = unit_resolution_binning_parameters(available_energies_including_zero(sc;negative_energies),sc.dims,sc;kwargs...)
+Base.copy(p::BinningParameters) = BinningParameters(copy(p.binstart),copy(p.binend),copy(p.binwidth),copy(p.covectors))
 
-function unit_resolution_binning_parameters(ωvals::AbstractVector{Float64})
-    if !all(abs.(diff(diff(ωvals))) .< 1e-12)
-      @warn "Non-uniform bins will be re-spaced into uniform bins"
-    end
-    if length(ωvals) == 1
-      error("Can not infer bin width given only one bin center")
-    end
-    ωbinwidth = (maximum(ωvals) - minimum(ωvals)) / (length(ωvals) - 1)
-    ωstart = minimum(ωvals) - ωbinwidth / 2
-    ωend = maximum(ωvals)
+# Support numbins as a (virtual) property, even though only the binwidth is stored
+Base.getproperty(params::BinningParameters, sym::Symbol) = sym == :numbins ? [count_bins(params.binstart[i],params.binend[i],params.binwidth[i]) for i = 1:4] : getfield(params,sym)
 
-    return ωstart, ωend, ωbinwidth
+function Base.setproperty!(params::BinningParameters, sym::Symbol, numbins)
+    if sym == :numbins
+        # *Ensure* that the last bin contains params.binend
+        params.binwidth .= (params.binend .- params.binstart) ./ (numbins .- 0.5)
+    else
+        setfield!(params,sym,numbins)
+    end
 end
 
 """
-    slice_2D_binning_parameter(sc::SampledCorrelations, cut_from_q, cut_to_q, cut_bins::Int64, cut_width::Float64; plane_normal = [0,0,1],cut_height = cutwidth)
-
-Creates [`BinningParameters`](@ref) which make a cut along one dimension of Q-space.
- 
-The x-axis of the resulting histogram consists of `cut_bins`-many bins ranging
-from `cut_from_q` to `cut_to_q`. 
-The width of the bins in the transverse direciton is controlled by `cut_width` and `cut_height`.
-
-The binning in the transverse directions is defined in the following way, which sets their normalization and orthogonality properties:
-
-    cut_covector = normalize(cut_to_q - cut_from_q)
-    transverse_covector = normalize(plane_normal × cut_covector)
-    cotransverse_covector = normalize(transverse_covector × cut_covector)
-
-In other words, the axes are orthonormal with respect to the Euclidean metric.
-
-If the cut is too narrow, there will be very few scattering vectors per bin, or
-the number per bin will vary substantially along the cut.
-If the output appears under-resolved, try increasing `cut_width`.
-
-The four axes of the resulting histogram are:
-  1. Along the cut
-  2. Fist transverse Q direction
-  3. Second transverse Q direction
-  4. Energy
-
-This function can be used without reference to a [`SampledCorrelations`](@ref) using this alternate syntax to manually specify the bin centers for the energy axis:
+    count_bins(binstart,binend,binwidth)
 
-    slice_2D_binning_parameter(ω_bincenters, cut_from, cut_to,...)
+Returns the number of bins in the binning scheme implied by `binstart`, `binend`, and `binwidth`.
+To count the bins in a [`BinningParameters`](@ref), use `params.numbins`.
 
-where `ω_bincenters` specifies the energy axis, and both `cut_from` and `cut_to` are arbitrary covectors, in any units.
+This function defines how partial bins are handled, so it should be used preferentially over
+computing the number of bins manually.
 """
-function slice_2D_binning_parameters(ωvals::Vector{Float64},cut_from_q,cut_to_q,cut_bins::Int64,cut_width;plane_normal = [0,0,1],cut_height = cut_width)
-    # This covector should measure progress along the cut in r.l.u.
-    cut_covector = normalize(cut_to_q - cut_from_q)
-    # These two covectors should be perpendicular to the cut, and to each other
-    transverse_covector = normalize(plane_normal × cut_covector)
-    cotransverse_covector = normalize(transverse_covector × cut_covector)
-
-    start_x = cut_covector ⋅ cut_from_q
-    end_x = cut_covector ⋅ cut_to_q
-
-    transverse_center = transverse_covector ⋅ cut_from_q # Equal to using cut_to_q
-    cotransverse_center = cotransverse_covector ⋅ cut_from_q
-
-    ωstart, ωend, ωbinwidth = unit_resolution_binning_parameters(ωvals)
-    xstart, xend, xbinwidth = unit_resolution_binning_parameters(range(start_x,end_x,length = cut_bins))
-
-    binstart = [xstart,transverse_center - cut_width/2,cotransverse_center - cut_height/2,ωstart]
-    binend = [xend,transverse_center,cotransverse_center,ωend]
-    numbins = [cut_bins,1,1,length(ωvals)]
-    covectors = [cut_covector... 0; transverse_covector... 0; cotransverse_covector... 0; 0 0 0 1]
-
-    BinningParameters(binstart,binend;numbins = numbins, covectors = covectors)
-end
-
-function slice_2D_binning_parameters(sc::SampledCorrelations,cut_from_q,cut_to_q,args...;kwargs...)
-    slice_2D_binning_parameters(available_energies_including_zero(sc),cut_from_q,cut_to_q,args...;kwargs...)
+function count_bins(bin_start,bin_end,bin_width)
+  if !isfinite(bin_width)
+    1
+  else
+    ceil(Int64,(bin_end - bin_start) / bin_width)
+  end
 end
 
 """
@@ -302,65 +258,213 @@ The following alternative syntax can be used to compute bin centers for a single
     axes_bincenters(binstart,binend,binwidth)
 """
 function axes_bincenters(binstart,binend,binwidth)
-    bincenters = Vector{AbstractRange{Float64}}(undef,0)
+    bincenters = Vector{AbstractVector{Float64}}(undef,0)
     for k = eachindex(binstart)
-        first_center = binstart[k] .+ binwidth[k] ./ 2
-        nbin = count_bins(binstart[k],binend[k],binwidth[k])
-        push!(bincenters,range(first_center,step = binwidth[k],length = nbin))
+        if isfinite(binwidth[k])
+          first_center = binstart[k] .+ binwidth[k] ./ 2
+          nbin = count_bins(binstart[k],binend[k],binwidth[k])
+          push!(bincenters,range(first_center,step = binwidth[k],length = nbin))
+        else
+          push!(bincenters,[binstart[k]])
+        end
     end
     bincenters
 end
 axes_bincenters(params::BinningParameters) = axes_bincenters(params.binstart,params.binend,params.binwidth)
 
 function axes_binedges(binstart,binend,binwidth)
-    binedges = Vector{AbstractRange{Float64}}(undef,0)
+    binedges = Vector{AbstractVector{Float64}}(undef,0)
     for k = eachindex(binstart)
-        nbin = count_bins(binstart[k],binend[k],binwidth[k])
-        push!(binedges,range(binstart[k],step = binwidth[k],length = nbin + 1))
+        if isfinite(binwidth[k])
+          nbin = count_bins(binstart[k],binend[k],binwidth[k])
+          push!(binedges,range(binstart[k],step = binwidth[k],length = nbin + 1))
+        else
+          push!(binedges,[-Inf,Inf])
+        end
     end
     binedges
 end
 axes_binedges(params::BinningParameters) = axes_binedges(params.binstart,params.binend,params.binwidth)
 
+# Find an axis-aligned bounding box containing the histogram
+function binning_parameters_aabb(params)
+    (; binstart, binend, covectors) = params
+    bin_edges = axes_binedges(params)
+    first_edges = map(x -> x[1],bin_edges)
+    last_edges = map(x -> x[end],bin_edges)
+    bin_edges = [first_edges last_edges]
+    this_corner = MVector{4,Float64}(undef)
+    q_corners = MMatrix{4,16,Float64}(undef)
+    for j = 1:16 # The sixteen corners of a 4-cube
+        for k = 1:4 # The four axes
+            this_corner[k] = bin_edges[k,1 + (j >> (k-1) & 1)]
+        end
+        this_corner[.!isfinite.(this_corner)] .= 0
+        q_corners[:,j] = covectors \ this_corner
+    end
+    lower_aabb_q = minimum(q_corners,dims=2)[1:3]
+    upper_aabb_q = maximum(q_corners,dims=2)[1:3]
+    return lower_aabb_q, upper_aabb_q
+end
 
+struct BinnedIntensities{T} <: AbstractIntensities
+    # Original chemical cell
+    #crystal :: Crystal
+    # BinningParameters in RLU
+    params :: BinningParameters
+    # Intensity data as bin-integrated values
+    data :: Array{T, 4} # collect(size(data)) == params.numbins
+    # Number of individually binned contributions (useful for some normalizations)
+    counts :: Array{Float64, 4}
+end
 
-"""
-    reciprocal_space_path_bins(sc,qs,density,args...;kwargs...)
+function Base.show(io::IO, res::BinnedIntensities)
+    sz = join(res.params.numbins, "×")
+    print(io, string(typeof(res)) * " ($sz bins)")
+end
 
-Takes a list of wave vectors, `qs` in R.L.U., and builds a series of histogram [`BinningParameters`](@ref)
-whose first axis traces a path through the provided points.
-The second and third axes are integrated over according to the `args` and `kwargs`,
-which are passed through to [`slice_2D_binning_parameters`](@ref).
+function binned_intensities(sc,params::BinningParameters;kT = nothing,integrated_kernel = nothing)
+    static_mode = sc isa SampledCorrelationsStatic
+    if !isnothing(integrated_kernel) && static_mode
+      error("Can't broaden if data is not energy-resolved")
+    end
 
-Also returned is a list of marker indices corresponding to the input points, and
-a list of ranges giving the indices of each histogram `x`-axis within a concatenated histogram.
-The `density` parameter is given in samples per reciprocal lattice unit (R.L.U.).
+    # Decide on which Q points can possibly contribute (depends on geometry of
+    # supercell and params)
+    lower_aabb_q, upper_aabb_q = binning_parameters_aabb(params)
+    # Round the axis-aligned bounding box *outwards* to lattice sites
+    # SQTODO: are these bounds optimal?
+    Ls = size((static_mode ? sc.parent : sc).data)[4:6]
+    lower_aabb_cell = floor.(Int64,lower_aabb_q .* Ls .+ 1) 
+    upper_aabb_cell = ceil.(Int64,upper_aabb_q .* Ls .+ 1) 
+    cells = CartesianIndices(Tuple(((:).(lower_aabb_cell,upper_aabb_cell))))[:]
+    qpts = QPoints([Vec3((cell.I .- 1) ./ Ls) for cell = cells])
+
+    # Calculate intensity at prepared qpts. No broadening yet because
+    # that depends on the bin layout and uses an integrated_kernel!
+    energies = if static_mode
+            [0.0]
+        else
+            sort(available_energies_including_zero(sc;negative_energies=true))
+        end
+    res = if static_mode
+      intensities_static(sc, qpts)
+    else
+      intensities(sc, qpts; energies, kT)
+    end
+
+    # Bin (and broaden) those intensities according to BinningParameters
+    k = MVector{3,Float64}(undef)
+    v = MVector{4,Float64}(undef)
+    q = view(v,1:3)
+    coords = MVector{4,Float64}(undef)
+    xyztBin = MVector{4,Int64}(undef)
+    xyzBin = view(xyztBin,1:3)
+
+    (; binwidth, binstart, binend, covectors, numbins) = params
+    return_type = typeof(res).parameters[1]
+    output_intensities = zeros(return_type,numbins...)
+    output_counts = zeros(Float64,numbins...)
+
+    # Pre-compute discrete broadening kernel from continuous one provided
+    if !isnothing(integrated_kernel)
+        # Upgrade to 2-argument kernel if needed
+        integrated_kernel_edep = try
+            integrated_kernel(0.,0.)
+            integrated_kernel
+        catch MethodError
+            (ω,Δω) -> integrated_kernel(Δω)
+        end
+
+        fraction_in_bin = Vector{Vector{Float64}}(undef,length(energies))
+        for (iω,ω) in enumerate(energies)
+            fraction_in_bin[iω] = Vector{Float64}(undef,numbins[4])
+            for iωother = 1:numbins[4]
+                ci_other = CartesianIndex(xyzBin[1],xyzBin[2],xyzBin[3],iωother)
+                # Start and end points of the target bin
+                a = binstart[4] + (iωother - 1) * binwidth[4]
+                b = binstart[4] + iωother * binwidth[4]
+
+                # P(ω picked up in bin [a,b]) = ∫ₐᵇ Kernel(ω' - ω) dω'
+                fraction_in_bin[iω][iωother] = integrated_kernel_edep(ω,b - ω) - integrated_kernel_edep(ω,a - ω)
+            end
+        end
+    end
+
+    for cell_ix = 1:length(cells), (iω,ω) in enumerate(energies)
+        cell = cells[cell_ix]
+        q .= ((cell.I .- 1) ./ Ls) # q is in R.L.U.
+        k .= (static_mode ? sc.parent : sc).crystal.recipvecs * q
+        if isnothing(integrated_kernel) # `Delta-function energy' logic
+            # Figure out which bin this goes in
+            v[4] = ω
+            mul!(coords,covectors,v)
+            coords .= (coords .- binstart) ./ binwidth
+            coords[.!isfinite.(binwidth)] .= 0
+            xyztBin .= 1 .+ floor.(Int64,coords)
+
+            # Check this bin is within the 4D histogram bounds
+            if all(xyztBin .<= numbins) && all(xyztBin .>= 1)
+                intensity = static_mode ? res.data[cell_ix] : res.data[iω,cell_ix]
+
+                ci = CartesianIndex(xyztBin.data)
+                output_intensities[ci] += intensity
+                output_counts[ci] += 1
+            end
+        else # `Energy broadening into bins' logic
+            # For now, only support broadening for `simple' energy axes
+            if covectors[4,:] == [0,0,0,1] && norm(covectors[1:3,:] * [0,0,0,1]) == 0
+
+                # Check this bin is within the *spatial* 3D histogram bounds
+                # If we are energy-broadening, then scattering vectors outside the histogram
+                # in the energy direction need to be considered
+                mul!(view(coords,1:3),view(covectors,1:3,1:3), view(v,1:3))
+                xyzBin .= 1 .+ floor.(Int64,(view(coords,1:3) .- view(binstart,1:3)) ./ view(binwidth,1:3))
+                if all(xyzBin .<= view(numbins,1:3)) &&  all(xyzBin .>= 1)
+
+                    # Calculate source scattering vector intensity only once
+                    intensity = res.data[iω,cell_ix]
+
+                    # Broaden from the source scattering vector (k,ω) to
+                    # each target bin ci_other
+                    ci_other = CartesianIndex(xyzBin[1],xyzBin[2],xyzBin[3])
+                    view(output_intensities,ci_other,:) .+= fraction_in_bin[iω] .* Ref(intensity)
+                    view(output_counts,ci_other,:) .+= fraction_in_bin[iω]
+                end
+            else
+                error("Energy broadening not yet implemented for histograms with complicated energy axes")
+            end
+        end
+    end
+    N_bins_in_BZ = abs(det(covectors[1:3,1:3])) / prod(binwidth[1:3])
+    output_data = output_intensities ./ N_bins_in_BZ ./ length(energies)
+    BinnedIntensities(params,output_data,output_counts)
+end
+
+"""
+    integrate_axes!(params::BinningParameters; axes)
+Integrate over one or more axes of the histogram by setting the number of bins
+in that axis to 1. Examples:
+
+    integrate_axes!(params; axes = [2,3])
+    integrate_axes!(params; axes = 2)
 """
-function q_space_path_bins(ωvals,qs,density,args...;kwargs...)
-    nPts = length(qs)
-    params = []
-    markers = []
-    ranges = []
-    total_bins_so_far = 0
-    push!(markers, total_bins_so_far+1)
-    for k = 1:(nPts-1)
-        startPt = qs[k]
-        endPt = qs[k+1]
-        # Density is taken in R.L.U. since that's where the
-        # scattering vectors are equally spaced!
-        nBins = round(Int64,density * norm(endPt - startPt))
-        # SQTODO: Automatic density that adjusts itself lower
-        # if there are not enough (e.g. zero) counts in some bins
-
-        param = slice_2D_binning_parameters(ωvals,startPt,endPt,nBins,args...;kwargs...)
-        push!(params,param)
-        push!(ranges, total_bins_so_far .+ (1:nBins))
-        total_bins_so_far = total_bins_so_far + nBins
-        push!(markers, total_bins_so_far+1)
+function integrate_axes!(params::BinningParameters;axes)
+    for k in axes
+        nbins = [params.numbins.data...]
+        nbins[k] = 1
+        params.numbins = SVector{4}(nbins)
     end
-    return params, markers, ranges
+    return params
+end
+
+function energy_resolve!(params::BinningParameters,energies)
+  energies = sort(energies)
+  params.binend[4] = maximum(energies)
+  params.binwidth[4] = energies[2] - energies[1]
+  params.binstart[4] = energies[1] - params.binwidth[4]/2
+  params
 end
-q_space_path_bins(sc::SampledCorrelations, qs::Vector, density,args...;kwargs...) = q_space_path_bins(available_energies_including_zero(sc), qs, density,args...;kwargs...)
 
 function available_energies_including_zero(x; kwargs...)
     ωs = available_energies(x;kwargs...)
@@ -368,4 +472,4 @@ function available_energies_including_zero(x; kwargs...)
     (length(ωs) == 1 && isnan(ωs[1])) ? [0.] : ωs
 end
 
-=#
+