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 x = 1 :  w - 1
      postep = floor((x)/(w)*100)
      for y = 2 : h - 1
          old = out(x, y, :);
          old(1,1,:) = max(squeeze(old), [0.0 0.0 0.0]');
          old(1,1,:) = min(squeeze(old), [255.0 255.0 255.0]');

          col = old(:)';
          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(x:x+1, y-1:y+1, :) = out(x:x+1, y-1:y+1, :) + errwin;
          out(x, y, :) = new;
      end
  end

 out = uint8(out);