This approach has several drawbacks. Not only does it quickly become very slow for high spatial resolutions (unless a shader is implemented), it also does not return the decision boundaries. Ideally, one would for each location $\mathbf{r}_k$ obtain a sorted list of vertices defining the boundary of the points belonging to class $k$. In fact, obtaining the decision boundaries in a voronoi diagram is equivalent to Delaunay triangulation -- finding triangles for a set of points such that the circumcirle of each triangle does not contain points from other triangles.
An interesting algorithm for obtaining the Voronoi boundaries for a set of locations is Fortunes algorithm. Armed with this algorithm one may compute the Voronoi boundaries linearly, $O(N)$, in the number of locations $N$. The idea is to use a line, referred to as a sweepline, and for each location, consider the curves that are equidistant to the location and the line -- parabolas. Let us, without loss of generality, assume this line to be vertical. Then the parabolas may be stiched together to form a piecewise parabolic curve by, for each horisontal line, considering only the parabola closest to the sweepline from the left side. This line, and the subset of locations involved, is referred to as the beach. A location is added to the beach when it is passed by the sweepline. A location is removed (from that specific position within the beach) when three consecutive parabolas in the beach intersect. At such an intersection point, being a point where three locations are equidistant to each other (and to the sweepline), are vertices in the Voronoi diagram.
When pressing (mouse-press) the canvas, you can see the fully rendered voronoi diagram. There are some bugs in the code, but most of the time it works. Since this is my playground, I did not feel like debugging this. Maybe I will at some point...