Lloyd's algorithm in normed vector spaces How do I run Lloyd's algorithm in a normed vector space? 
The space:
L*a*b* color space, finite sRGB segment, $R^3$
The distance metric:
CIE94 using L*C*h* information derived from the L*a*b* coordinates.
The sub-problems (so far):


*

*how to compute a voronoi diagram in a normed vector space.

*how to compute the centroid and therefore volume in such a space.


If I am working with a polygon (ie, a voronoi cell) whose points are computed relative to the previous point using the distance metric, I don't believe I can represent the points contained within the polygon using straight lines. In other words, while the triangle inequality may hold (single shortest path between two points), the shortest path may not be a straight line. This means, as far as I know, I can't depend on line-based algorithms to calculate voronoi diagrams. Furthermore this complicates calculating the volume of geometric primitives such as spheres and prisms.
Additionally, do the concepts of volume (necessary for calculating the centre of mass) carry-over from Euclidean space? For example, is the volume of the rectangular prism formed from two points on opposing corners, p1 and p2, still calculated as: d(p1.x, p2.x) * d(p1.y - p2.y) * d(p1.z - p2.z).
Alternatively:
I am working only with 2**24 different points - one point for each unique colour in a 3-channel 8-bit unsigned integer image. When working with voxels (a unit cube, one per point) in Euclidean space the voronoi diagram can be created in O(n) time (n == # of voxels) using a flood-fill algorithm, which has benefit of dynamically updating the centroid for each voronoi cell (the centre of mass is updated each time a voxel is added)
As an alternative to working in a non-Euclidean space, I would accept any answer that demonstrates how to map points in the CIE94 space back to the L*a*b* (Euclidean) space preserving both distance and volume. That is, for each of the 2**24 points p in L*a*b* space, map p to a new point p1 in L*a*b* space that preserves volume and distance per the CIE94 distance metric.
Overview:
My goal is to find the N most perceptually distinguishable colours, given P pre-existing colours. Since the colour difference in L*a*b* is equivalent to distance, this problem generalizes to the sphere-packing problem in an irregular (sRGB) space.
To handle P pre-existing colours I will modify Lloyd's algorithm such that only the points in N will be shifted.
I think (probably) that Lloyd's algorithm will only get me a Random Close Pack configuration. Therefore the starting positions of the points in N will be as closely aligned to the centres of spheres in a Hexagonal Close Pack arrangement as the points in P will allow. However, this will only be possible using the alternative solution (mapping); I can't do this in the CIE94-norm vector space.
Update: Optimizations
As has been mentioned below, LLoyd's algorithm can be generalized over a set of points using only the distance metric. While Lloyd's algorithm inherently generates a voronoi diagram, in my question I've failed to discern between the algorithm itself, and algorithms for computing voronoi diagrams in order to optimize the speed of Lloyd's algorithm.
These are some algorithms I have considered:
Voronoi + Centroid: Calculate the voronoi diagram via Delauny triangulation, and the centroid by integrating the points forming the boundaries of the voronoi cells. Complexity is $O(M*log(M)) + O(M) == O(nlogn)$, where M represents the number of voronoi cells, or reproduction points in Lloyd's algorithm. This method works in any infinite space (for example the bounded $R^3$ sRGB space), so I can extend my algorithm to handle an arbitrary number of unique colours, rather than simply 2**24. This would be useful if I ever have to process floating-point images. Unfortunately, I don't know how to extend this to non-Euclidean spaces.
Floodfill algorithm: Calculate the voronoi diagram by flooding the voxel grid (N voxels, one per unique colour, in my case 2**24) outwards from each point M (the number of voronoi cells, or reproduction points in Lloyd's algorithm. Complexity is $O(n)$. however this method is only applicable to a set of points and not generic spaces. This can't be extended to non-Euclidean spaces (without varying the size of a voxel).
Therefore my question is more appropriately phrased as: How do I efficiently run Lloyd's algorithm in $R^3$ with on either a set of points or a bounded space, when the distance metric is not Euclidean?
 A: You do not need remapping or any such step provided that you know your distance (distortion) function and have a way to generate samples of your original distribution. Lloyd algorithm is a pretty powerful algorithm that can be applied (in practice) to any distance function that one can imagine. For example, consider an arbitrary distance function $d(x,\hat{x})$, where $x$ represents your original sample, and $\hat{x}$ is its reproduction. $x$ can be a scalar, vector, matrix, whatever you want (your color coordinates), and in fact $x$ and $\hat{x}$ may even lie in different spaces. All of these are irrelevant for the discussion below, which is general. 
Now, suppose your original samples have a PDF $f(x)$. In practice, you generate a sufficiently large amount of training vectors, say $x_1,\ldots,x_N$, where $N$ is large. At this stage, you also decide how precise you want to be in the process of reconstructing the $x$s, and suppose you will allow $M$ reproduction points $\hat{x}_1,\ldots,\hat{x}_M$. Then, the Lloyd algorithm will work as follows:
1) Assign each $x_i$ to its corresponding reproduction point that will minimize the distance, and call this assignment, say $\alpha_i\in\{1,\ldots,M\}$. We have
$\alpha_i = \arg\min_{j\in\{1,\ldots,M\}} d(x_i,\hat{x}_j)$.
2) Recalculate the reproduction points $\hat{x}_i = \arg\min_{z} \sum_{j:\alpha_j = i}d(x_j, z)$, and loop.
Step 1 is trivial and will not pose any difficulty. For step 2, if we were working with squared-error distance, say $d(x,\hat{x}) = \|x-\hat{x}\|^2$, then we would have $\hat{x}_i = \frac{\sum_{j:\alpha_j = i} x_j}{|\{j:\alpha_j = i\}|}$. For a general $d$, you will not have such a neat closed-form solution to the $\arg\min$ problem in the second step, and therefore you will have to resort to some optimization tools which will give you the optimal $z$, but this is also a minor headache in my opinion. 
