function out = fsd(in, pal)
  [h, w, z] = size(in);
  p_h = size(pal, 1);

  out = zeros(h, w, z, "double");
  out(1:h-1, 2:w-1, :) = double(in(1:h-1, 2:w-1, :));

  pal = double(pal);

  win = [0 0 7; 3 5 1]/16;

  for j = 1 : h - 1
      for i = 2 :  w - 1
          old = out(j, i, :);

          col = old(:)';
          col = max(col, [0.0 0.0 0.0]);
          col = min(col, [255.0 255.0 255.0]);
          col = repmat(col, p_h, 1);

          tmp = pal - col;
          tmp = tmp .^ 2;
          tmp = sum(tmp, 2);

          [c, idx] = min(tmp);
          new = cat(3, pal(idx, 1), pal(idx, 2), pal(idx, 3));

          err = old - new;
          errwin = cat(3, err(:,:,1) * win, err(:,:,2) * win, err(:,:,3) * win);

          out(j : j + 1, i - 1 : i + 1, :) = out(j : j + 1, i - 1 : i + 1, :) + errwin;
          out(j, i, :) = new;
      end
  end

 out = uint8(out);