Yes, the algorithm is correct. It can be sped up by using implementing the 'count number of points within a box' using smarter data structures. In particular, a running time of $O(n^2 \log n)$ is achievable.
Let $N(x_\ell,y_\ell,x_u,y_u)$ denote the number of points within the box $[x_\ell,x_u] \times [y_\ell,y_u]$, and define $N'(x_u,y_u) := N(-\infty,-\infty,x_u,y_u)$. Note that
$$N(x_\ell,y_\ell,x_u,y_u) = N'(x_u,y_u) - N'(x_u,y_\ell) - N'(x_\ell,y_u) + N'(x_\ell,y_\ell).$$
I will describe below two ways to compute $N'(\cdot,\cdot)$ efficiently, which immediately leads to more efficient versions of your algorithm.
Store all points in a 2D range tree, 2D interval tree, or 2D segment tree. I'm not sure, but I think these data structures allow you to efficiently compute $N'(\cdot,\cdot)$, though I'm not sure if there are any worse-case running time bounds.
If you care about provable worse-case complexity, here is an approach that will enable you to solve your problem in $O(n^2 \log n)$ time, using a sweep line algorithm and persistent data structures (a standard technique for algorithms on 2D points).
Build a self-balancing binary tree with one leaf per unique $y$-coordinate of the points (thus at most $n$ leaves). We sweep a vertical line from left to right. We'll have one tree per $x$-coordinate of the points (with repetition). Each node of the tree corresponds to a range of $y$-coordinates (correspondings to the leaves that are descendants of that node). We'll augment the tree for $x$-coordinate $x_u$ so that the node corresponding to range $[y_\ell,y_u]$ stores the count $N(-\infty,y_\ell,x_u,y_u)$.
Start with the leftmost $x$-coordinate, and build the tree for this $x$-coordinate. Then, as you sweep the line to the right by one position, you pick up one additional point $(x_u,y_u)$, so you need to increment the count for the leaf corresponding to $(x_u,y_u)$, as well as ancestors of that leaf. Instead of modifying the tree in place, use a persistent data structure to implement this modification, so we retain the tree before this modification and the tree afterwards; these trees differ in about $\lg n$ nodes and the remaining $n-\lg n$ nodes are common to both trees. Continue sweeping the sweep line to the right, modifying the tree as you encounter each point.
Once you have built all the trees, you can compute $N'(x_u,y_u)$ by looking up the tree for $x$-coordinate $x_u$, then finding the set of internal nodes whose ranges form a partition of $[-\infty,y_u]$. This a set of $O(\log n)$ nodes. Summing up the counts in those nodes gives us the value of $N'(x_u,y_u)$.
What is the running time? It takes $O(n \log n)$ time to build the persistent data structure with all the trees. Also, to compute $N'(x_u,y_u)$, you can find the tree for a particular point, find the partition of nodes, and sum up all of them in $O(\log n)$ time. Therefore, computing $N(x_\ell,y_\ell,x_u,y_u)$ can be done in $O(\log n)$.
Therefore, the total running time of the algorithm becomes $O(n^2 \log n)$, as you have to compute $N(\cdot)$ for $n^2$ pairs of points.
I suspect perhaps there is a clever algorithm to achieve $O(n^2)$ running time, but I can't see how at the moment.