Make half-size spanning trees in the RGB cube and the pixels, then form hamiltonian cycles through both. This general idea has been done before, but no source code exists online so I had to figure out the relevant algorithms myself. What hasn't been done before is biasing the weights used for the spanning trees, which has the effect of forming interesting patterns on the final image while still keeping the hamiltonian property. Takes about 10s to generate a 4096x4096 image. Code repo:
https://github.com/mnurzia/hamrgb