I wanted to achieve the best color and shading reproduction with a single-pass algorithm. This image is converted in about 30 seconds using NumPy and SciKit-Image, and I know that further optimizations are possible.
The source image and a full-RGB image are converted to YIQ color space. The Q (purple-green) channel is not very useful, so it is replaced by another Y (brightness) channel to accentuate shading. Each pixel is binned into a channel and sorted, then replaced with the matching pixel from the full RGB image.
