Linear gradient (C)

From LiteratePrograms
(Redirected from Linear gradient)
Jump to: navigation, search

This article will describe a way to efficiently draw a linear gradient, using control points similar to a drawing program: given two points and a list of colors (we will refer them as color-stop afterward), draw a series of lines perpendicular to the initial segment, interpolating colors along the path.

The mathematical background behind this algorithm are digital differential analyzer (DDA) and a little bit of trigonometry. The implementation will be written in C99 and using SDL as a rendering backend.

Contents

[edit] First approach

Given the constraints in the introduction, we can already define the prototype of our function, along a few datatypes:

<<gradient_v1.h>>=

#ifndef LINEAR_GRADIENT_H
#define LINEAR_GRADIENT_H

Image data structure

ColorStop data structure

int LinearGradient(Image img, int x1, int y1, int x2, int y2, ColorStop cs, int count);
#endif

The Image data structure stores the dimensions of the image and the stride (bytes per scanline), which may differ from width*bpp/8 due to alignment (this kind of datatype is also useful if you want to only alter a smaller rectangular portion of the bitmap):

<<Image data structure>>=
typedef struct Image_t *       Image;
typedef uint8_t *              DATA8;

struct Image_t
{
	int   width, height, bpp;
	int   stride;
	DATA8 bitmap;
};


ColorStop will describe which colors the gradient is made of. Obviously, at least two colors are required. dist is an optional parameter to describe at which position the color begins. We will use a number between 0 and 100, that we will scale accordingly at the initialization stage.

<<ColorStop data structure>>=
typedef struct ColorStop_t *   ColorStop;

struct ColorStop_t
{
	uint8_t rgba[4];
	int     dist;
};

The approach taken to render the gradient will be by drawing it by scanline/column (depending on the orientation of the initial segment). Drawing the gradient by rendering directly the perpendicular lines using whatever line rendering algorithm (like Bresenham) usually produce very poor results because of the staircase effect (the algorithm we will describe will also be subject to this effect, but it will be much easier to correct).

[edit] Trigonometry

Figure 1: input parameters.

To render the gradient by scanline, we will need to project each color stop on the first scanline. This can be easily done using a bit of trigonometry. To better visualize what's going on, check figure on the right.

There is two parameters we need to find the value of: x0 and g. x0 can be easily found using some trigonometry within the right triangle ADx0 (angles are in degrees):

\tan(90-\alpha) = \frac{y_1}{AD}

AD = \frac{y_1}{\tan(90-\alpha)}

\tan(90-\alpha) = \mathrm{cotan}\ \alpha = \frac{x_2-x_1}{y_2-y_1}

x_0 = x_1 + AD = x_1 + y_1 \frac{y_2-y_1}{x_2-x_1}

Finding g, we can use some properties of the right triangle ABC:

\cos\ \alpha = \frac{AB\ \mathrm{(adjacent)}}{g\ \mathrm{(hypotenuse)}}

g = \frac{AB}{cos\ \alpha}

AB can be easily found using the following formula:

\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}

The cosine of the α angle can be more easily found by considering the triangle ABE:

\cos\ \alpha = \frac{AE}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}

Therefore the value of g will be :

g = AB \frac{AB}{AE} = \frac{(x_2-x_1)^2+(y_2-y_1)^2}{x_2-x_1}

Which can computed by:

<<compute g>>=
adj = x2 - x1;
opp = y2 - y1;
g = (adj*adj+opp*opp) / adj;

Knowing those results we can already see that x2 = x1 is a special case (vertical gradient). Actually to simplify even further the explanations, we will divide the algorithm into octant, according to the angle made by segment AB in respect to the horizontal axis. In a first approach, we will only see what to expect in the first octant, that is between the angle 0 and 45 degrees (remember that the origin is the top left corner, as it is usual for computer screen). Later we will see that all the other octants can be rendered like the first by slightly adjusting input parameters.

Knowing x0 and g, we can now easily render the gradient: we have to iterate through the scanline of the image, interpolating the corresponding values made by the line (0,x0) − (x1,y1), then draw a gradient from the first color stop to the last color stop over the next g pixels. The pixels before x0 will be set to the same color than the first color stop and the pixels after x0 + g will be set to the same color than the last color stop.

[edit] Digital Differential Analyzer

Digital Differential Analyzer (DDA) can, given two points (x1, y1) and (x2, y2), iterate over one axis and linearly interpolate corresponding values on the other axis, using very simple operation (mostly integer addition and subtraction). They can be seen as a special case of the Bresenham algorithm, but instead of giving all the integer points there are in a segment, DDA will only give you one point per axis you iterate on. The math behind them is pretty straightforward: given our two points, we can deduce the amount the y coordinate increase each time the x coordinate is increased:

slope = Δy / Δx

This is a floating number, not the ideal type when summed iteratively. So instead of using a floating number, we separate the division into quotient and remainder; the quotient is the (integer) amount the y coordinate increase for each increase of the x coordinate, and the remainder is the error accumulated. The error range goes from 0 (including) to Δx (not including). Whenever the error accumulated exceed Δx, we have to increase the y coordinate again.

Let's wrap this into some helper functions:

<<Iter structure>>=
typedef struct Iter_t          Iter;

struct Iter_t
{
	int x, y, xe, ye;
	int dx, dy, err, sx, sy, oldy;
};

<<DDA helper functions>>=
void InitDDA(Iter * iter, int xs, int xe, int ys, int ye)
{
	/* Pre-condition: xe > xs >= 0 */
	div_t q = div(iter->dy = ye - ys, xe);
	iter->y   = ys;
	iter->ye  = ye;
	iter->x   = xs;
	iter->xe  = xe;
	iter->err = xe >> 1;
	iter->dx  = abs(q.rem);
	iter->sx  = q.quot;
	iter->sy  = (ys < ye ? 1 : -1);
	iter->oldy = -1;
	if (xs > 0)
	{
		q = div(xs * iter->dy + (xe >> 1), xe);
		iter->y   = ys + q.quot;
		iter->err = xe - q.rem;
	}
}

static inline void IterDDA(Iter * iter)
{
	iter->oldy = iter->y;
	iter->x ++;
	iter->y += iter->sx;
	iter->err -= iter->dx;
	if (iter->err < 0) iter->y += iter->sy, iter->err += iter->xe;
}
#define ISDDAEND(iter) ((iter).x >= (iter).xe)

Some note here :

  • The error accumulated is initially set to Δx/2, to do automatic up or down rounding.
  • The axis we are supposed to iterate is x, starting at (0, ys), ending at (xe, ye). xs is to start at another point on this line (useful for clipping). The code path for starting at another point is not very optimized, because there will be optimizations that will make this case extremely unlikely: for drawing a gradient this code path might be executed at most once.
  • Instead of adding the error, we subtract it. This slightly simplify the test when resetting the error.
  • Nothing prevents you from iterating past xe. We will use this trick when enumerating coordinates from the perpendicular.

[edit] First version

Knowing all this, we can implement the first version of our linear gradient drawing routine.

The main LinearGradient routine loops over the scanlines, drawing each one with DrawScanline, using a DDA to track the starting column of each one:

<<main gradient loop>>=
Iter pos;
InitDDA(&pos, 0, y1, x1 + y1 * opp / adj, x1);
while (pos.x < img->height)
{
	DrawScanline(img, pos.y, pos.x, cs, count, g);
	IterDDA(&pos);
}

The rest of LinearGradient, preceding this loop, does parameter validation and scales the cs.dist values to the range [0,g]:

<<LinearGradient initialization>>=
int adj, opp, g, i;

/* Check for sane parameters */
if (x1 == x2 && y1 == y2) return 0;
if (count < 2) return 0;

compute g

/* Adjust range of cs.dist for each color from [0, 100] to [0, g] */
for (i = 0; i < count; i ++)
{
	int dist;
	if (i == 0) dist = 0; else
	if (i == count-1) dist = g; else
	if (cs[i].dist == 0) dist = i * g / (count-1);
	else dist = cs[i].dist * g / 100;
	cs[i].dist = dist;
}

We can then define LinearGradient as:

<<LinearGradient definition>>=
int LinearGradient(Image img, int x1, int y1, int x2, int y2, ColorStop cs, int count)
{
	LinearGradient initialization

	main gradient loop

	return 1;
}

The DrawScanline function uses a separate DDA to track the current color of each color channel, skipping any pixels outside the frame:

<<DrawScanline declaration>>=
static void DrawScanline(Image img, int x, int y, ColorStop cs, int count, int max);

<<SETU32 helper macro>>=
#define SETU32(dst, src)      * (uint32_t *) (dst) = * (uint32_t *) (src)

<<DrawScanline definition>>=
/* Draw scanline 'y', starting gradient at column 'x', over 'max' pixels */
static void DrawScanline(Image img, int x, int y, ColorStop cs, int count, int max)
{
	DATA8 p = img->bitmap + img->stride * y;
	ColorStop c = cs;
	Iter r, g, b, a;
	int i, w = img->width;

	if (x > 0) {
		/* Set first iter->x px to color c0 */
		uint8_t rgb[4];
		memcpy(rgb, c->rgba, 4);
		for (i = x; i > 0; SETU32(p, rgb), p += 4, i --);
	}

	/* Skip first -x px */
	if (x < 0) {
		for (c ++, i = 1; -x > c->dist && i < count; c ++, i ++);
		if (i < count) x += c[-1].dist;
		i = -x; x = 0;
	}
	else i = 0, c++;

	/* Draw gradient */
	do {
		int len = c->dist - c[-1].dist;
		InitDDA(&r, i, len, c[-1].rgba[3], c->rgba[3]);
		InitDDA(&g, i, len, c[-1].rgba[2], c->rgba[2]);
		InitDDA(&b, i, len, c[-1].rgba[1], c->rgba[1]);
		InitDDA(&a, i, len, c[-1].rgba[0], c->rgba[0]);

		do {
			p[3] = r.y; IterDDA(&r);
			p[2] = g.y; IterDDA(&g);
			p[1] = b.y; IterDDA(&b);
			p[0] = a.y; IterDDA(&a); x ++; p += img->bpp;
			if (ISDDAEND(r)) {
				if (c->dist == max) {
					/* Keep using the last color, up to the end of line */
					r.sx = r.err = r.dx = 0;
					g.sx = g.err = g.dx = 0;
					b.sx = b.err = b.dx = 0;
					a.sx = a.err = a.dx = 0;
					r.xe = 1<<30;
					if (x < w) continue;
				}
				else c ++;
				break;
			}
		}
		while (x < w);
	} while (x < w);
}

The main() function, the program entry point, initializes SDL, calls LinearGradient, and runs an event loop:

<<main entry point>>=
int main (int argc, char * argv[])
{
	if (SDL_Init(SDL_INIT_VIDEO) < 0)
		return 1;

	atexit(SDL_Quit);

	SDL_Surface * screen = SDL_SetVideoMode(500, 400, 32, SDL_HWSURFACE|SDL_ANYFORMAT);
	if (! screen) return 1;

	struct Image_t img = {
		.bitmap = screen->pixels, .stride = screen->pitch, .width = screen->w, .height = screen->h,
		.bpp = screen->format->BytesPerPixel
	};

	/* Simple red to blue gradient */
	#define  NBCOL     2
	struct ColorStop_t colors[NBCOL] = {{{255, 0, 0}}, {{0, 20, 255}}};
	int i;

	/* We need to convert colors into native format */
	for (i = 0; i < NBCOL; i ++)
	{
		DATA8    rgba = colors[i].rgba;
		uint32_t col = SDL_MapRGBA(screen->format, rgba[0], rgba[1], rgba[2], rgba[3]);
		SETU32(colors[i].rgba, &col);
	}

	LinearGradient(&img, 50, 100, 450, 350, colors, NBCOL);
	SDL_UpdateRect(screen, 0, 0, screen->w, screen->h);

	while (1)
	{
		SDL_Event event;
		while (SDL_WaitEvent(&event))
		{
			switch (event.type) {
			case SDL_QUIT: return 0;
			case SDL_KEYDOWN:
				if (event.key.keysym.sym == SDLK_ESCAPE) return 0;
			}
		}
	}
	return 0;
}

Finally we put it all together to get our first working gradient program:

<<gradient_v1.c>>=
#include <stdlib.h>
#include <SDL/SDL.h>
#include "gradient_v1.h"

Iter structure

DDA helper functions

DrawScanline declaration

LinearGradient definition

SETU32 helper macro

DrawScanline definition

main entry point

Image generated by the above code.

In Linux this program can be built with gcc, using sdl-config to determine the right compiler flags:

<<build_and_run_v1.sh>>=
gcc gradient_v1.c -o gradient_v1 `sdl-config --cflags --libs`
./gradient_v1

The result should be the image to the right.

[edit] Horizontal predictor

Before we see how to generalize this algorithm to other octants, there is an optimization that can drastically reduce the number of operations. When we iterate over the perpendicular line, we notice that each line is shifted by at most 1 pixel, compared to the previous line. This is because the rendering is limited to one octant, thus the perpendicular angle varies from -90 to -45 degrees. That means, after drawing the first scanline, each subsequent ones can be drawn by copying the previous one, shifting it by 1, 0 or -1 pixel, and possibly iterating the gradient where the previous scanline stopped. So, in the worst case scenario, we will have to compute one scanline and one column, all the other pixels will be copied from previous lines, eliminating almost all of DDA iteration.

However, this optimization will modify quite a bit of the original algorithm and have some side effects when we will generalize it later, but the speedup is well worth the effort.

There is two optimizations we'll do : keep the iterating state across scanlines and use the horizontal predictor to skip most of computing operations. To do this, we need to introduce a new datatype that will keep that drawing state:

<<GradState data structure>>=
typedef struct GradState_t     GradState;

struct GradState_t
{
	Iter r, g, b, a, pos;
	int step, dir, err;
	ColorStop c;
};

When we introduced DDA, there was a field that did not have an obvious usage: oldy. This is where it will become handy.

<<main gradient loop with GradState>>=
/* Draw the gradient: init DDA to scan perpendicular values */
GradState state;
state.c = cs;
state.step = count;
InitDDA(&state.pos, 0, y1, x1 + y1 * opp / adj, x1);
while (state.pos.x < img->height)
{
	DrawScanline(img, &state, g);
	IterDDA(&state.pos);
}

The rest of LinearGradient is unmodified:

<<LinearGradient with GradState definition>>=
int LinearGradient(Image img, int x1, int y1, int x2, int y2, ColorStop cs, int count)
{
	LinearGradient initialization

	main gradient loop with GradState

	return 1;
}

DrawScanline is modified to use GradState as well:

<<DrawScanline with GradState declaration>>=
static void DrawScanline(Image img, GradState * state, int max);

<<DrawScanline with GradState definition>>=
static void DrawScanline(Image img, GradState * state, int max)
{
	DATA8 p = img->bitmap + img->stride * state->pos.x;
	ColorStop c = state->c;
	int i, diff, w = img->width, chan = img->bpp;
	int x = state->pos.y;

	/* Horizontal predictor */
	if (state->pos.x > 0)
	{
		int size;
		diff = state->pos.y - state->pos.oldy;
		w    = abs(diff);
		size = (img->width - w) * chan;
		if (diff < 0) {
			memcpy(p, p - img->stride + chan, size);
			p += (img->width - 1) * chan;
		}
		else memcpy(p + chan * diff, p - img->stride, size);
		/* No differences between this line and the previous one */
		if (w == 0) return;
		x = 0;
	}
	else /* First line */
	{
		if (x > 0)
			/* Set first iter->x px to color c0 */
			for (i = x; i > 0; SETU32(p, c->rgba), p += 4, i --);

		/* Skip first -x px */
		if (x < 0) {
			for (c ++, i = state->step-1; -x > c->dist && i > 0; c ++, i --);
			x += c[-1].dist;
			i = -x; x = 0;
		}
		else i = 0, c++;

		diff = c->dist - c[-1].dist;
		InitDDA(&state->r, i, diff, c[-1].rgba[3], c->rgba[3]);
		InitDDA(&state->g, i, diff, c[-1].rgba[2], c->rgba[2]);
		InitDDA(&state->b, i, diff, c[-1].rgba[1], c->rgba[1]);
		InitDDA(&state->a, i, diff, c[-1].rgba[0], c->rgba[0]);
	}

	/* Draw gradient */
	do {
		p[3] = state->r.y; IterDDA(&state->r);
		p[2] = state->g.y; IterDDA(&state->g);
		p[1] = state->b.y; IterDDA(&state->b);
		p[0] = state->a.y; IterDDA(&state->a); x ++; p += 4;
		if (ISDDAEND(state->r)) {
			if (c->dist == max) {
				/* Keep using the last color, up to the end of line */
				state->r.sx = state->r.err = state->r.dx = 0;
				state->g.sx = state->g.err = state->g.dx = 0;
				state->b.sx = state->b.err = state->b.dx = 0;
				state->a.sx = state->a.err = state->a.dx = 0;
				state->r.xe = 1<<30;
				if (x < w) continue; else break;
			}
			else c ++;
			/* Next color stop */
			diff = c->dist - c[-1].dist;
			InitDDA(&state->r, 0, diff, c[-1].rgba[3], c->rgba[3]);
			InitDDA(&state->g, 0, diff, c[-1].rgba[2], c->rgba[2]);
			InitDDA(&state->b, 0, diff, c[-1].rgba[1], c->rgba[1]);
			InitDDA(&state->a, 0, diff, c[-1].rgba[0], c->rgba[0]);
		}
	} while (x < w);
	state->c = c;
}

This produces our second implementation, which is otherwise identical to the first:

<<gradient_with_GradState.c>>=
#include <stdlib.h>
#include <SDL/SDL.h>
#include "gradient_v1.h"

Iter structure

DDA helper functions

GradState data structure

DrawScanline with GradState declaration

LinearGradient with GradState definition

SETU32 helper macro

DrawScanline with GradState definition

main entry point

<<build_and_run_with_GradState.sh>>=
gcc gradient_with_GradState.c -o gradient_with_GradState `sdl-config --cflags --libs`
./gradient_with_GradState

The output is perfectly identical to the first version's.

[edit] Anti-aliasing

Staircase effect for nearly vertical lines. Click to see in detail.

One last minor modification we will see before generalizing this algorithm, is how to avoid the "staircase effect" (see right). This happens because when enumerating the values of the perpendicular line, we are rounding the values to the nearest integer. Since this line is almost vertical, there will be a whole lot of scanlines with the same iter.y value, creating this "staircase effect".

To avoid this, we can use the accumulated error within the DDA, which can be seen as a number between 0 and 1 (iter.err / iter.xe). This is actually the amount the scanline need to be horizontally shifted to match its real position on the screen. Since this amount is less than one pixel, we will have to do sub-pixel anti-aliasing. Also one minor adjustment, is that the accumulated error is initialized to Δx/2. Since we will offset the error, we do not need that automatic up/down rounding (although, we still want it for iterating over the gradient values), therefore the accumulated error has to be initialized to Δx.

In this section, we will just briefly explain the strategies, along a few code samples, because the algorithm will have to be significantly modified again for the generalized version.

First, we will define a function that will shift a scanline, according to a normalized value (16-bit unsigned int).

<<AntiAlias: performs subpixel antialiasing>>=
static void AntiAlias(Image img, int err, int pos)
{
	DATA8 p = img->bitmap + pos;
	int   c = img->bpp;
	int   i = img->width-1;
	if (err == 0) return;
	for (p += i*c, c = -c; i > 0; i --, p += c) {
		p[0] = ((65536 - err) * p[0] + err * p[c]   + 65536/2) >> 16;
		p[1] = ((65536 - err) * p[1] + err * p[c+1] + 65536/2) >> 16;
		p[2] = ((65536 - err) * p[2] + err * p[c+2] + 65536/2) >> 16;
		p[3] = ((65536 - err) * p[3] + err * p[c+3] + 65536/2) >> 16;
	}
}
Result of using different rounding methods with anti-aliasing and rendering the gradient. Click to see in more detail.

When applying anti-aliasing, you have to be careful to use the same rounding method than the one used to render the gradient, otherwise you will get artifacts, as shown right (faint horizontal and periodic break). The horizontal white lines on the left have been added to spot on which scanlines the artifacts are. You might have to adjust the gamma of your monitor to see them. Also the image has been magnified 4 times.

Next step is to normalize the error from the DDA. We have to distinguish positive from negative slope (state.pos.sy). If the slope is positive, the error has to be added to current interpolated value. If it is negative, it has to be substracted and we need to subtract 1 from interpolated value, to avoid doing left/right shift in the anti-aliasing phase.

Here is also where the horizontal predictor gets in the way: we cannot apply anti-aliasing to the line we just draw, because the horizontal predictor will need it for the next scanline. Therefore, we will have to slightly delay the anti-aliasing, by using a loop like this (excerpt from LinearGradient()):

/* Draw the gradient: init DDA to scan perpendicular values */
GradState state;
int err[2], pos = 0, error;
state.c = cs;
state.step = count;
InitDDA(&state.pos, 0, y1, x1 + y1 * opp / adj, x1);
state.pos.err = state.pos.xe;

while (state.pos.x < img->height)
{
	/* Normalize the error from the DDA */
	if (state.pos.sy < 0) error = (state.pos.err << 16) / state.pos.xe;
	else                  error = ((state.pos.xe - state.pos.err) << 16) / state.pos.xe;

	err[state.pos.x&1] = error;

	DrawScanline(img, &state, g);
	if (state.pos.x > 0)
		AntiAlias(img, error = err[(state.pos.x-1)&1], pos, state.pos.sy), pos += img->stride;
	IterDDA(&state.pos);
}
return 1;

[edit] Generalization

To identify octant, we will use two parameters: Δx = x2-x1 and Δy = y2-y1.

[edit] Horizontal segment (|Δx| > |Δy|)

Δx>0 and Δy>0
This is the octant we study so far.
Δx<0 and Δy<0
Exchange x2 and x1, y2 and y1, and reverse color stops of the gradient. Then render the gradient like the first octant.
Δx>0 and Δy<0
Mirror x2 and x1, and render the gradient from right to left (do not reverse color stops).
Δx<0 and Δy>0
Mirror x2 and x1, then draw the gradient like octant 2, but reversed (right to left).

[edit] Vertical segment (|Δx| < |Δy|)

There are 2 ways to handle vertical octants: either we duplicate the DrawScanline() function to draw pixel by column instead of by line. Or, we keep the original DrawScanline() and rotate the scanline on the fly. Since this function is not that trivial, we will do the latter. So the overall strategy here is to allocate a temporary image of 2 scanlines (one for rendering the gradient and one for the horizontal predictor) of img->height pixels, and then write a custom anti-aliasing function that will put those pixels back into the column of the original image, while applying anti-aliasing. The function AntiAliasV that does this is similar to the AntiAlias function seen earlier:

<<AntiAliasV for vertical segments>>=
/* Copy scanline into column and apply anti-aliasing */
static void AntiAliasV(Image img, int err, int pos, DATA8 src)
{
	DATA8 p = img->bitmap + pos;
	int   c = img->stride, inc = - img->bpp;
	int   i = img->height-1;
	for (p += i*c, c = -c, src += i * -inc; i > 0; i --, p += c, src += inc) {
		p[0] = ((65536 - err) * src[0] + err * src[inc]   + 65536/2) >> 16;
		p[1] = ((65536 - err) * src[1] + err * src[inc+1] + 65536/2) >> 16;
		p[2] = ((65536 - err) * src[2] + err * src[inc+2] + 65536/2) >> 16;
		p[3] = ((65536 - err) * src[3] + err * src[inc+3] + 65536/2) >> 16;
	}
	SETU32(p, src);
}

[edit] Final version

Knowing all this, we can put together the final version of the gradient function, working in all 8 octants:

<<gradient_general.h>>=
#ifndef LINEAR_GRADIENT_H
#define LINEAR_GRADIENT_H

Image data structure

ColorStop data structure

Iter structure

GradState data structure

int LinearGradient(Image img, int x1, int y1, int x2, int y2, ColorStop cs, int count);

#ifndef SETU32
SETU32 helper macro
#endif
#ifndef swap
#define swap(a, b) do { int temp = a; a = b; b = temp; } while(0)
#endif

#endif

This chunk is large or complex. Large chunks are challenging for readers to follow. Consider breaking it up into more pieces with more prose explanation.

<<gradient_general.c>>=
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <malloc.h>
#include "gradient_general.h"

DDA helper functions

static void DrawScanline(Image img, GradState * state, int max);
static void AntiAlias(Image img, int err, int pos);
static void AntiAliasV(Image img, int err, int pos, DATA8 src);

int LinearGradient(Image img, int x1, int y1, int x2, int y2, ColorStop cs, int count)
{
	int adj, opp, g, i, horz, inc;
	ColorStop colors, c;
	GradState state;

	/* Check for sane parameters */
	if (x1 == x2 && y1 == y2) return 0;
	if (count < 2) return 0;

	adj = x2 - x1;
	opp = y2 - y1;
	horz = abs(adj) >= abs(opp);
	state.dir = img->bpp;
	if (! horz) { swap(opp, adj); swap(x1, y1); swap(x2, y2); }
	if (opp * adj < 0) {
		state.dir = - state.dir;
		i = horz ? img->width : img->height;
		x1 = i - x1; x2 = i - x2;
		adj = x2 - x1;
	}
	g = (adj*adj+opp*opp) / adj;
	c = colors = alloca(sizeof *c * count);
	inc = 1;
	state.c = c;
	state.step = count;
	if (g < 0) c += count-1, inc = -1, c->dist = abs(g);
	else colors->dist = 0;

	/* Adjust range of cs.dist for each color from [0, 100] to [0, g] */
	for (i = 0; i < count; i ++, c += inc)
	{
		int dist;
		SETU32(c->rgba, cs[i].rgba);
		if (i == 0) dist = 0; else
		if (i == count-1) dist = g; else
		if (cs[i].dist == 0) dist = i * g / (count-1);
		else dist = cs[i].dist * g / 100;
		if (inc > 0)
		{
			if (dist < 0) dist = -dist;
		} else {
			if (g < 0) dist -= g;
			else dist = g - dist;
		}
		c->dist = dist;
	}

	/* Draw the gradient: init DDA to scan perpendicular values */
	if (y1 == y2)
	{
		/* Horizontal gradient, handle special case where y1 == y2 == 0 */
		InitDDA(&state.pos, 0, 1, x1, x1);
	}
	/* Use ref point the farther away from y = 0, for extra precision */
	else if (abs(y1) < abs(y2))
	{
		int x0 = x1 + y1 * opp / adj;
		InitDDA(&state.pos, 0, y2, x0, x0 - opp * y2 / adj);
	}
	else InitDDA(&state.pos, 0, y1, x1 + y1 * opp / adj, x1);
	state.pos.err = state.pos.xe;
	state.err = state.dir * state.pos.sy;
	if (horz)
	{
		int err[2], pos = 0, error;
		while (state.pos.x < img->height)
		{
			if (state.err < 0) error = (state.pos.err << 16) / state.pos.xe;
			else               error = ((state.pos.xe - state.pos.err) << 16) / state.pos.xe;
			err[state.pos.x&1] = error;

			DrawScanline(img, &state, g);
			if (state.pos.x > 0)
				AntiAlias(img, error = err[(state.pos.x-1)&1], pos), pos += img->stride;
			IterDDA(&state.pos);
		}
	}
	else /* iterate over vertical axis */
	{
		int   chan = img->bpp, max = img->width, error, pos;
		Image imgtmp = malloc(sizeof *img + img->height * chan * 2);
		DATA8 line;
		if (imgtmp == NULL) return 0;

		memcpy(imgtmp, img, sizeof *img);
		imgtmp->bitmap = line = (DATA8) (imgtmp+1);
		imgtmp->width  = img->height;
		imgtmp->height = img->width;
		imgtmp->stride = img->height * chan;
		pos = 0;
		do {
			if (state.err < 0) error = (state.pos.err << 16) / state.pos.xe;
			else               error = ((state.pos.xe - state.pos.err) << 16) / state.pos.xe;
			DrawScanline(imgtmp, &state, g);
			AntiAliasV(img, error, pos, line); pos += chan;
			if (state.pos.x > 0)
				memcpy(imgtmp->bitmap, line = imgtmp->bitmap + imgtmp->stride, imgtmp->stride);
			IterDDA(&state.pos); state.pos.x = 1; max --;
		}
		while (max > 0);
		free(imgtmp);
	}
	return 1;
}

AntiAlias: performs subpixel antialiasing

AntiAliasV for vertical segments

/* Draw scanline 'pos.x', starting gradient at column 'pos.y', over 'max' pixels */
static void DrawScanline(Image img, GradState * state, int max)
{
	DATA8 p = img->bitmap + img->stride * state->pos.x;
	ColorStop c = state->c;
	int i, diff, w = img->width, chan = img->bpp;
	int x = state->pos.y;

	if (max < 0) x += max, max = -max;
	if (state->pos.sy < 0 && state->pos.err > 0) x --;

	/* Horizontal predictor */
	if (state->pos.x > 0)
	{
		int size;
		diff = state->pos.y - state->pos.oldy;
		w    = abs(diff);
		size = (img->width - w) * chan;
		if (state->dir < 0) diff = - diff;
		if (diff < 0) {
			memcpy(p, p - img->stride + chan, size);
			p += (img->width - 1) * chan;
		}
		else memcpy(p + chan * diff, p - img->stride, size);
		/* No differences between this line and the previous one */
		if (w == 0) return;
		x = 0; chan = state->dir;
	}
	else /* First line */
	{
		chan = state->dir;
		if (chan < 0) p += (img->width - 1) * -chan;
		if (x > 0)
			/* Set first iter->x px to color c0 */
			for (i = x; i > 0; SETU32(p, c->rgba), p += chan, i --);

		/* Skip first -x px */
		if (x < 0) {
			for (c ++, i = state->step-1; -x > c->dist && i > 0; c ++, i --);
			x += c[-1].dist;
			i = -x; x = 0;
		}
		else i = 0, c++;

		diff = c->dist - c[-1].dist;
		InitDDA(&state->r, i, diff, c[-1].rgba[3], c->rgba[3]);
		InitDDA(&state->g, i, diff, c[-1].rgba[2], c->rgba[2]);
		InitDDA(&state->b, i, diff, c[-1].rgba[1], c->rgba[1]);
		InitDDA(&state->a, i, diff, c[-1].rgba[0], c->rgba[0]);
	}

	/* Draw gradient */
	do {
		p[3] = state->r.y; IterDDA(&state->r);
		p[2] = state->g.y; IterDDA(&state->g);
		p[1] = state->b.y; IterDDA(&state->b);
		p[0] = state->a.y; IterDDA(&state->a); x ++; p += chan;
		if (ISDDAEND(state->r)) {
			if (c->dist == max) {
				/* Keep using the last color, up to the end of line */
				state->r.sx = state->r.err = state->r.dx = 0;
				state->g.sx = state->g.err = state->g.dx = 0;
				state->b.sx = state->b.err = state->b.dx = 0;
				state->a.sx = state->a.err = state->a.dx = 0;
				state->r.xe = 1<<30;
				if (x < w) continue; else break;
			}
			else c ++;
			/* Next color stop */
			diff = c->dist - c[-1].dist;
			InitDDA(&state->r, 0, diff, c[-1].rgba[3], c->rgba[3]);
			InitDDA(&state->g, 0, diff, c[-1].rgba[2], c->rgba[2]);
			InitDDA(&state->b, 0, diff, c[-1].rgba[1], c->rgba[1]);
			InitDDA(&state->a, 0, diff, c[-1].rgba[0], c->rgba[0]);
		}
	} while (x < w);
	state->c = c;
}

This chunk is large or complex. Large chunks are challenging for readers to follow. Consider breaking it up into more pieces with more prose explanation.

<<gradtest.c>>=
#include <SDL/SDL.h>
#include <math.h>
#include "gradient_general.h"

#define WINW   500
#define WINH   400
#define NBCOL  2
static struct ColorStop_t colors[NBCOL] = {{{255, 0, 0}}, {{20, 20, 255}}};
SDL_Surface * screen;
int angle;

static Uint32 push_ticks(Uint32 interval, void *param)
{
	SDL_Event event = {.type = SDL_USEREVENT};
	SDL_PushEvent(&event);
	return interval;
}

static void reset(Image img)
{
	angle ++;

	int x = cos(angle * M_PI / 180) * 220;
	int y = sin(angle * M_PI / 180) * 220;
	memset(img->bitmap, 0, img->stride * img->height);
	LinearGradient(img, WINW/2-x, WINH/2-y, WINW/2+x, WINH/2+y, colors, NBCOL);
	SDL_UpdateRect(screen, 0, 0, screen->w, screen->h);
}

int main (int argc, char * argv[])
{
	if (SDL_Init(SDL_INIT_VIDEO|SDL_INIT_TIMER) < 0)
		return 1;

	atexit(SDL_Quit);

	screen = SDL_SetVideoMode(WINW, WINH, 32, SDL_HWSURFACE|SDL_ANYFORMAT);
	if (! screen) return 1;

	struct Image_t img = {
		.bitmap = screen->pixels, .stride = screen->pitch, .width = screen->w, .height = screen->h,
		.bpp = screen->format->BytesPerPixel
	};

	int i;

	/* We need to convert colors into native format */
	for (i = 0; i < NBCOL; i ++)
	{
		DATA8    rgba = colors[i].rgba;
		uint32_t col = SDL_MapRGBA(screen->format, rgba[0], rgba[1], rgba[2], rgba[3]);
		SETU32(colors[i].rgba, &col);
	}

	SDL_AddTimer(10, push_ticks, NULL);

	while (1)
	{
		SDL_Event event;
		while (SDL_WaitEvent(&event))
		{
			switch (event.type) {
			case SDL_USEREVENT: reset(&img); break;
			case SDL_QUIT: return 0;
			case SDL_KEYDOWN:
				if (event.key.keysym.sym == SDLK_ESCAPE) return 0;
			}
		}
	}
	return 0;
}

#ifdef WIN32
#include <windows.h>
int WINAPI WinMain (HINSTANCE instance,
                    HINSTANCE previnst,
                    LPSTR args,
                    int wndState)
{
	return main(0, NULL);
}
#endif

<<build_and_run_general.sh>>=
gcc gradtest.c gradient_general.c -o gradient_general `sdl-config --cflags --libs`
./gradient_general

Download code
hijacker
hijacker
hijacker
hijacker