I order all of the colors by walking a 3-dimensional Hilbert Curve via an L-System. I then walk the pixels in the output image along a 2-dimensional Hilbert Curve and lay out all of the colors.

512 x 512 output:

Here's the Scala code. Most of it covers just the logic and math of moving through three dimensions via pitch/roll/yaw. I'm sure there was a better way to do that part, but oh well.

import scala.annotation.tailrec
import java.awt.image.BufferedImage
import javax.imageio.ImageIO
import java.io.File

object AllColors {
  
  case class Vector(val x: Int, val y: Int, val z: Int) {
    def applyTransformation(m: Matrix): Vector = {
      Vector(m.r1.x * x + m.r1.y * y + m.r1.z * z, m.r2.x * x + m.r2.y * y + m.r2.z * z, m.r3.x * x + m.r3.y * y + m.r3.z * z)
    }
    def +(v: Vector): Vector = {
      Vector(x + v.x, y + v.y, z + v.z)
    }
    def unary_-(): Vector = Vector(-x, -y, -z)
  }
  
  case class Heading(d: Vector, s: Vector) {
    def roll(positive: Boolean): Heading = {
      val (axis, b) = getAxis(d)
      Heading(d, s.applyTransformation(rotationAbout(axis, !(positive ^ b))))
    }
    
    def yaw(positive: Boolean): Heading = {
      val (axis, b) = getAxis(s)
      Heading(d.applyTransformation(rotationAbout(axis, positive ^ b)), s)
    }
    
    def pitch(positive: Boolean): Heading = {
      if (positive) {
        Heading(s, -d)
      } else {
        Heading(-s, d)
      }
    }
    
    def applyCommand(c: Char): Heading = c match {
      case '+' => yaw(true)
      case '-' => yaw(false)
      case '^' => pitch(true)
      case 'v' => pitch(false)
      case '>' => roll(true)
      case '<' => roll(false)
    }
  }
  
  def getAxis(v: Vector): (Char, Boolean) = v match {
    case Vector(1, 0, 0) => ('x', true)
    case Vector(-1, 0, 0) => ('x', false)
    case Vector(0, 1, 0) => ('y', true)
    case Vector(0, -1, 0) => ('y', false)
    case Vector(0, 0, 1) => ('z', true)
    case Vector(0, 0, -1) => ('z', false)
  }
  
  def rotationAbout(axis: Char, positive: Boolean) = (axis, positive) match {
    case ('x', true) => XP
    case ('x', false) => XN
    case ('y', true) => YP
    case ('y', false) => YN
    case ('z', true) => ZP
    case ('z', false) => ZN
  }
  
  case class Matrix(val r1: Vector, val r2: Vector, val r3: Vector)
  
  val ZP = Matrix(Vector(0,-1,0),Vector(1,0,0),Vector(0,0,1))
  val ZN = Matrix(Vector(0,1,0),Vector(-1,0,0),Vector(0,0,1))
  
  val XP = Matrix(Vector(1,0,0),Vector(0,0,-1),Vector(0,1,0))
  val XN = Matrix(Vector(1,0,0),Vector(0,0,1),Vector(0,-1,0))
  
  val YP = Matrix(Vector(0,0,1),Vector(0,1,0),Vector(-1,0,0))
  val YN = Matrix(Vector(0,0,-1),Vector(0,1,0),Vector(1,0,0))
  
  @tailrec def applyLSystem(current: Stream[Char], rules: Map[Char, List[Char]], iterations: Int): Stream[Char] = {
    if (iterations == 0) {
      current
    } else {
      val nextStep = current flatMap { c => rules.getOrElse(c, List(c)) }
      applyLSystem(nextStep, rules, iterations - 1)
    }
  }
  
  def walk(x: Vector, h: Heading, steps: Stream[Char]): Stream[Vector] = steps match {
    case Stream() => Stream(x)
    case 'f' #:: rest => x #:: walk(x + h.d, h, rest)
    case c #:: rest => walk(x, h.applyCommand(c), rest)
  }
  
  def hilbert3d(n: Int): Stream[Vector] = {
    val rules = Map('x' -> "^>x<f+>>x<<f>>x<<+fvxfxvf+>>x<<f>>x<<+f>x<^".toList)
    val steps = applyLSystem(Stream('x'), rules, n) filterNot (_ == 'x')
    walk(Vector(0, 0, 0), Heading(Vector(1, 0, 0), Vector(0, 1, 0)), steps)
  }
  
  def hilbert2d(n: Int): Stream[Vector] = {
    val rules = Map('a' -> "-bf+afa+fb-".toList, 'b' -> "+af-bfb-fa+".toList)
    val steps = applyLSystem(Stream('a'), rules, n) filterNot (c => c == 'a' || c == 'b')
    walk(Vector(0, 0, 0), Heading(Vector(1, 0, 0), Vector(0, 0, 1)), steps)
  }
  
  def main(args: Array[String]): Unit = {
    val n = 4
    val img = new BufferedImage(1 << (3 * n), 1 << (3 * n), BufferedImage.TYPE_INT_RGB)
    hilbert3d(n * 2).zip(hilbert2d(n * 3)) foreach { case (Vector(r,g,b), Vector(x,y,_)) => img.setRGB(x, y, (r << (24 - 2 * n)) | (g << (16 - 2 * n)) | (b << (8 - 2 * n))) }
    ImageIO.write(img, "png", new File(s"out_$n.png"))
  }
}