Efficiently generating 3D histograms for large sets 2D coordinates

Question

I'm attempting to generate a 3D histogram for a large set of two-dimensional coordinates ($\approx 10^7$). the x- and y-components of each coordinate are rounded to the nearest integer, and these rounded values are used as indices in a histogram matrix counting the number of coordinates rounded to the same integer values.

For example, if we have some list of elements:

ElementList={{102.2134124213,101.2421321312},{500.2,2761.7},{102.542423,101.344}}

We would add $+2$ to position {102, 101} in the histogram matrix, and $+1$ to the position {500, 2761}.

Here's the piece of code I'm running:

ElementList = Table[{RandomReal[{1, 2000}], RandomReal[{1, 2000}]}, {x, 1, 10^7}];

HistogramDimX = 2000;
HistogramDimY = HistogramDimX;

ImageMatrix = Array[0 &, Length[ElementList]];
ImageMatrixHistogram = Array[0 &, {HistogramDimX, HistogramDimY}];

For[i = 1, i <= Length[ElementList], i++,

RoundedCoordinate = Round[ElementList[[i]]];
ImageMatrixHistogram[[Round[RoundedCoordinate[[1]]], 
Round[RoundedCoordinate[[2]]]]] += 1;

];

MatrixPlot[ImageMatrixHistogram, ColorFunctionScaling -> True, MaxPlotPoints -> 10^12]

This takes 94.7 seconds to run on a single X5690 3.47 GHz Intel Xeon(R) CPU. Is there a way to significantly speed this process up, and use a more efficient data structure that scales proportionally with the number of elements in ElementList rather than the dimensions of the histogram matrix (i.e. HistogramDimX & HistogramDimY)? I suppose I'd very much like the output to look something like the output of MatrixPlot in the above example.

Directly applying MatrixPlot to the data works terribly, which I suppose is due to some auto-interpolation occurring.