How to fix the dithering method to create a color gradient?

Question

There is a program that builds a matrix for a color gradient from white to black. Then, a dithering algorithm is applied to the matrix to eliminate the "stripes". 4 methods of dithering are implemented: Ordered, Random, Floyd-Steinberg, Jarvis-Judice-Ninke. First, I create a matrix of a certain size, convert it to a gradient, and output the result to a file format .pgm, type P5. If I translate the file into .png, I get the following image:

However, when you zoom in on the image, you can see the stripes (if you look closely):

This is the result of the program without dithering. The problem is that if you apply one of the dithering algorithms to the matrix, then stripes remain on the image. The result suggests that dithering does not work. What could be wrong? Do I need to use dithering first, and then build a gradient? Or the error is that you need to create a matrix of float or double type? How can I fix it?

Code:

#include "stdafx.h"
#include <iostream>
#include<algorithm>
#include<iterator>
#include<fstream>
#include<vector>
#include<cassert>
#include <ctime>
#include <sstream>

using namespace std;

vector<vector<int>> make_gradient(int height, int width)
{
    assert(height > 0 && width > 0);

    int cf = height / 255;
    int color = 0;
    vector<vector<int>> result(height, vector<int>(width));
    for (int i = 0; i < height; i += cf)
    {
        for (int j = 0; j < cf; ++j)
        {
            fill(result[i + j].begin(), result[i + j].end(), color % 255);
        }
        color++;
    }
    stable_sort(result.begin(), result.end());
    return result;
}

vector<vector<int>> ordered_dither(int height, int width, vector<vector<int>> result)
{
    int ditherSize = 4;
    int diterLookup[] = { 0, 8, 2, 10, 12, 4, 14, 6, 3, 11, 1, 9, 15, 7, 13, 5 };

    for (int i = 0; i < height; i++)
    {
        for (int j = 0; j < width; j++)
        {
            int xlocal = i%ditherSize;
            int ylocal = j%ditherSize;
            int requiredShade = diterLookup[xlocal + ylocal * 4] * 255 / 16;
            if ((requiredShade > 0) && (requiredShade < 1))
            {
                if (requiredShade >= (result[i][j] % 1)) {
                    result[i][j] = floor(result[i][j]);
                }
                else {
                    result[i][j] = ceil(result[i][j]);
                }
            }
            else requiredShade = 0;
        }
    }
    return result;
}

vector<vector<int>> random_dither(int height, int width, vector<vector<int>> result)
{
    for (int i = 0; i < height; i++)
    {
        for (int j = 0; j < width; j++)
        {
            int requiredShade = (float)rand() / RAND_MAX;
            if ((requiredShade > 0) && (requiredShade < 1))
            {
                if (requiredShade >= (result[i][j] % 1)) {
                    result[i][j] = floor(result[i][j]);
                }
                else {
                    result[i][j] = ceil(result[i][j]);
                }
            }
            else requiredShade = 0;
        }
    }
    return result;
}

vector<vector<int>> fs_dither(int height, int width, vector<vector<int>> result)
{
    for (int i = 0; i < height; i++)
    {
        for (int j = 0; j < width; j++)
        {
            int oldpixel = result[i][j];
            int newpixel = round(result[i][j]);
            result[i][j] = newpixel;
            int quanterror = oldpixel - newpixel;
            if (j < width - 1) {
                result[i][j + 1] += quanterror * 7 / 16;
            }
            if (i < height - 1) {
                if (j > 0) {
                    result[i + 1][j - 1] += quanterror * 3 / 16;
                }
                result[i + 1][j] += quanterror * 5 / 16;
                if (j < width - 1) {
                    result[i + 1][j + 1] += quanterror * 1 / 16;
                }
            }
        }
    }
    return result;
}

vector<vector<int>> jjn_dither(int height, int width, vector<vector<int>> result)
{
    for (int i = 0; i < height; i++)
    {
        for (int j = 0; j < width; j++)
        {
            int oldpixel = result[i][j];
            int newpixel = round(result[i][j]);;
            result[i][j] = newpixel;
            int quanterror = oldpixel - newpixel;
            if (j < width - 1) {
                result[i][j + 1] += quanterror * 7 / 48;
                if (j<width - 2)
                    result[i][j + 2] += quanterror * 5 / 48;
            }

            if (i < height - 1) {
                if (j > 0) {
                    if (j > 1)
                        result[i + 1][j - 2] += quanterror * 3 / 48;
                    result[i + 1][j - 1] += quanterror * 5 / 48;
                }

                result[i + 1][j] += quanterror * 7 / 48;
                if (j < width - 1) {
                    result[i + 1][j + 1] += quanterror * 5 / 48;
                    if (j < width - 2)
                        result[i + 1][j + 2] += quanterror * 3 / 48;
                }
            }

            if (i < height - 2) {
                if (j > 0) {
                    if (j>1)
                        result[i + 2][j - 2] += quanterror * 1 / 48;
                    result[i + 2][j - 1] += quanterror * 3 / 48;
                }
                result[i + 2][j] += quanterror * 5 / 48;
                if (j < width - 1) {
                    result[i + 2][j + 1] += quanterror * 3 / 48;
                    if (j < width - 2)
                        result[i + 2][j + 2] += quanterror * 1 / 48;
                }
            }

        }
    }
    return result;
}

int main(int argc, char *argv[])
{
    if (argc < 5) {
        cout << "usage:" << endl << "prog.exe <filename> <width> <height> <dithering>" << endl;
        return 0;
    }
    stringstream w(argv[2]);
    stringstream h(argv[3]);
    stringstream d(argv[4]);
    int numcols, numrows, dithering;

    if (!(w >> numcols)) {
        cout << "width is not a number" << endl;
        return 0;
    }
    if (numcols < 1) {
        cout << "width must be more than zero" << endl;
        return 0;
    }

    if (!(h >> numrows)) {
        cout << "height is not a number" << endl;
        return 0;
    }
    if (numrows < 1) {
        cout << "height must be more than zero" << endl;
        return 0;
    }

    if (!(d >> dithering)) {
        cout << "dithering is not a number" << endl;
        return 0;
    }
    if (dithering < 0 || dithering>4) {
        cout << "dithering must be [0-4]" << endl;
        return 0;
    }

    srand(time(0));
    ofstream file;

    file.open(argv[1]);

    if (!file)
    {
        cout << "can't open file" << endl;
        return 0;
    }

    file << "P5" << "\n";

    file << numrows << " " << numcols << "\n";

    file << 255 << "\n";

    vector<vector<int>> pixmap{ make_gradient(numrows, numcols) };
    switch (dithering) {
    case 1:
        pixmap = ordered_dither(numrows, numcols, pixmap);
        break;
    case 2:
        pixmap = random_dither(numrows, numcols, pixmap);
        break;
    case 3:
        pixmap = fs_dither(numrows, numcols, pixmap);
        break;
    case 4:
        pixmap = jjn_dither(numrows, numcols, pixmap);
        break;
    default:
        break;
    }
    for_each(pixmap.begin(), pixmap.end(), [&](const auto& v) {
        copy(v.begin(), v.end(), ostream_iterator<char>{file, ""});
    });

    file.close();

}

Answer 1

It's interesting to see you using dither to get rid of those bands that you can hardly see -- in the old days we would only dither when we had to render in 4 bits per channel or so.

Anyway... Your first problem is that before you use dithering to reduce your gradient to 256 levels, you have to render it at more than 256 levels. make_gradient should probably be rendering the gradient at 65536 levels or even floating point.

Your second problem is, it seems to me, that your dithering is currently doing nothing at all. result[i][j] is an integer, so when you say something like result[i][j] = floor(result[i][j]); (I guess you ignore the compiler warnings about the conversion), it does nothing. If you were generating your gradient in floats that problem would go away too.

If you fix these problems, then your dither can work, however none of these dither methods is really optimal for operating at such closely spaced levels. There may still be banding artifacts of some kind remaining when you're done (although you'd have to look really closely to see them). To make your results look as good as possible, you should really use a TPDF dither with amplitude equal to two quantization levels. With coarsely space levels, this looks noisier than some of your other choices, but is statistically more uniform and will look better when the levels are finely spaced.

It's also easy -- just add two random numbers between -0.5 and 0.5 to each pixel before quantizing to integers.

TPDF is mentioned here: https://en.wikipedia.org/wiki/Dither , but it's most important because it's the type of dither used in sampling for signal processing to make sure that the quantization doesn't cause any first or second order artifacts.

EDIT:

I appreciate that you've been working on this, so here's the code to create the dithered gradient in final form in one step:

vector<vector<int>> make_dithered_gradient(int height, int width)
{
    assert(height > 0 && width > 0);

    vector<vector<int>> result(height, vector<int>(width));
    for (int i = 0; i < height; ++i)
    {
        // the target shade for each line of pixels is the average
        // ideal gradient value in that line, which is the same as the
        // ideal value in the middle of the line
        double target = ((double)i+0.5)*255.0/height;
        for (int j = 0; j < width; ++j)
        {
            double dither = ((double)rand()-(double)rand())/RAND_MAX;
            int val = (int)round(target+dither);
            if (val < 0)
                val = 0;
            if (val > 255)
                val = 255;
            result[i][j] = val;
        }
    }
    return result;
}