I got better results just dithering the rgb channels separately (so effectively an 8 colour palette, black, white, rgb, yellow, cyan, magenta). In p5js:

    var img
    var pixel
    var threshold
    var error = [0, 0, 0]
    var a0

    function preload() {
      img = loadImage("https://upload.wikimedia.org/wikipedia/commons/thumb/4/44/Albrecht_D%C3%BCrer_-_Hare%2C_1502_-_Google_Art_Project.jpg/1920px-Albrecht_D%C3%BCrer_-_Hare%2C_1502_-_Google_Art_Project.jpg")
    }

    function setup() {
      // I'm just using a low discrepancy sequence for a quasirandom
      // dither and diffusing the error to the right, because it's
      // trivial to implement
      a0 = 1/sqrt(5)
      pixelDensity(2)
      createCanvas(400, 400);
      image(img, 0, 0, 400, 400)
      loadPixels()
      pixel = 0
      threshold = 0
    }

    function draw() {
      if (pixel > 400*400*16) {
        return
      }
      for (var i = 0; i < 2000; i++) {
        threshold = (threshold + a0)%1
        for(var j=0; j< 3; j++) {
          var c = pixels[pixel + j]
          pixels[pixel + j] = c + error[j] > threshold * 255 ? 255 : 0
          error[j] += c - pixels[pixel + j]
        }
        pixel += 4
      }
      updatePixels()
    }

Of course this isn't trying to pick the closest colour in the palette as you're doing - it's just trying to end up with the same intensity of rgb as the original image. It does make me wonder if you should be using the manhattan distance instead of euclidean, to get the errors to add correctly.